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THE a( 3) SCHEME— A FOURTH-ORDER SPACE-TIME FLUX-CONSERVING 
AND NEUTRALLY STABLE CESE SOLVER 

Sin-Chung Chang 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 

Abstract 

The CESE development is driven by a belief that a solver should (i) enforce conservation laws in both 
space and time, and (ii) be built from a non-dissipative (i.e., neutrally stable) core scheme so that the 
numerical dissipation can be controlled effectively. To provide a solid foundation for a systematic CESE 
development of high order schemes, in this paper we describe the a(3) scheme — a new 4th-order space-time 
flux-conserving and neutrally stable CESE solver of the advection equation du/dt + adu/dx = 0. The space- 
time stencil of this two-level explicit scheme is formed by one point at the upper time level and three points 
at the lower time level. Because it is associated with three independent mesh variables u" , (u x )j, and ( u xx )" 
(the numerical analogues of u, du/dx, and d 2 u/dx , respectively) and three equations per mesh point, the 
new scheme is referred to as the a(3) scheme. As in the case of other similar CESE neutrally stable solvers, 
the a(3) scheme enforces conservation laws in space-time locally and globally, and it has the basic, forward 
marching, and backward marching forms. These forms are equivalent and satisfy a space-time inversion 
invariant property which is shared by the advection equation. (In physics, space-time inversion invariance is 
referred to as PT invariance where P denotes a parity, i.e., mirror-image or spatial-reflection, operation while 
T denote a time-reversal operation.) Based on the concept of PT invariance, a set of algebraic relations 
involving the coefficient matrices of the a(3) scheme is developed. As it turns out, in the von Neumann 
analysis, these relations lead to the conclusion that the a(3) scheme must be neutrally stable when it is 
stable. Also, in the same analysis, it is proved rigorously that: (i) all three amplification factors (i.e., the 
eigenvalues of the 3x3 amplification matrix) of the a(3) scheme are of unit magnitude for all phase angles 
9 of the Fourier modes considered in the von Neumann analysis if and only if \v\ < 1/2 (v = a At /ax); (ii) 
the a(3) scheme is stable if and only if \v\ < 1/2; and (iii) the a(3) scheme is linearly unstable (in a sense 
to be defined) if \i/\ = 1/2. These theoretical results have been confirmed numerically. Moreover, through 
numerical experiments, it is established that the a(3) scheme generally is (i) 4th-order accurate for the mesh 
variables u" and {u x )™\ and (ii) 2nd-order accurate for (u xx )/j. However, in some exceptional cases, the 
scheme can achieve perfect accuracy aside from round-off errors. Finally the phase errors of the principal 
amplification factor of the a(3) scheme are evaluated numerically and shown to be 0(9 4 ), a sharp reduction 
from those of the a scheme (the original neutrally stable CESE solver) which are O(0 2 ). 
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1. Introduction 


The space-time conservation element and solution element (CESE) method is a high-resolution and 
genuinely multidimensional method for solving conservation laws [1-73]. Its nontraditional features include: 

(i) a unified treatment of space and time; (ii) the introduction of conservation elements (CEs) and solution 
elements (SEs) as the vehicles for enforcing space-time flux conservation; (iii) a novel time marching strategy 
that has a space-time staggered stencil at its core and, as such, fluxes at an interface can be evaluated 
without using any interpolation or extrapolation procedure (which, in turn, leads to the method’s ability 
to capture shocks without using Riemann solvers); (iv) the requirement that each scheme be built from a 
non-dissipative core scheme and, as a result, the numerical dissipation can be controlled effectively; and (v) 
the fact that mesh values of the physical dependent variables and their spatial derivatives are considered as 
independent marching variables to be solve for simultaneously. Note that CEs are non-overlapping space- 
time subdomains introduced such that (i) the computational domain can be filled by these subdomains; and 

(ii) flux conservation can be enforced over each of them and also over the union of any combination of them. 
On the other hand, SEs are space-time subdomains introduced such that (i) the boundary of each CE can 
be divided into several component parts with each of them belonging to a unique SE; and (ii) within a SE, 
any physical flux vector is approximated using simple smooth functions. In general, a CE does not coincide 
with a SE. 

Without using flux-splitting or other special techniques, since its inception in 1991 [1], the unstructured- 
mesh compatible CESE method has been used to obtain numerous accurate ID, 2D and 3D steady and 
unsteady flow solutions with Mach numbers ranging from 0.0028 to 10 [51] . The physical phenomena modeled 
include traveling and interacting shocks, acoustic waves, vortex shedding, viscous flows, detonation waves, 
cavitation, flows in fluid film bearings, heat conduction with melting and/or freezing, electrodynamics, MHD 
vortex, hydraulic jump, crystal growth, and chromatographic problems [3-73]. In particular, its unexpected 
simple non-reflecting boundary conditions [9,68] and rather unique capability to resolve both strong shocks 
and small disturbances (e.g., acoustic waves) simultaneously [13,15,16] makes the CESE method an effective 
tool for attacking computational aeroacoustics (CAA) problems. Note that the fact that the second-order 
CESE schemes can solve CAA problems accurately is an exception to the commonly-held belief that a 
second-order scheme is not adequate for solving CAA problems. Also note that, while numerical dissipation 
is needed for shock capturing, it may also result in annihilation of small disturbances. Thus a solver that can 
handle both strong shocks and small disturbances simultaneously must be able to overcome this difficulty. 

In spite of its nontraditional features and potent capabilities, the core ideas of the CESE method are 
simple. In fact, all of its key features are the inescapable results of an honest pursuit driven by these 
simple ideas. The first and foremost is the belief that the method must be solid in physics. As such, in 
the CESE development, conservation laws are enforced locally and globally in their natural space-time unity 
forms for ID, 2D and 3D cases. Moreover, because direct physical interaction generally occurs only among 
the immediate neighbors, use of the simplest stencil also becomes a CESE requirement. Obviously, this 
requirement is also very helpful in simplifying boundary-condition implementation. 

The second idea emerges from the realization that stability and accuracy are two competing issues in 
time-accurate computations, i.e. , too much numerical dissipation will degrade accuracy while too little of it 
will cause instability. In other words, to meet both accuracy and stability requirements, computation must be 
performed away from the edge (“cliff”) of instability but not too far from it. This represents a real dilemma 
in numerical method development. As an example, schemes with high-order accuracy generally have high 
accuracy and low numerical dissipation. However, they are susceptible to instability. In fact, in dealing with 
complicated real-world problems, stability of these schemes often is difficult to maintain without resorting to 
ad hoc treatments. To confront this issue head-on, in CESE development, generally it is required that a solver 
be built from a non-dissipative (i.e., neutrally stable) core scheme. By definition, computations involving 
a neutrally stable scheme are performed right on the edge of instability and therefore the numerical results 
generated are non-dissipative. As such numerical dissipation can be controlled effectively if the deviation of 
a solver from its non-dissipative core scheme can be adjusted using some built-in parameters. Note that the 
above idea also plays an essential role in the recent successful development of a family of Courant number 
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insensitive schemes [59,61,64,65,67]. 

Other CESE ideas are: (i) the flux at an interface be evaluated in a simple and consistent manner; 
(ii) genuinely multidimensional schemes be built as simple, consistent and straightforward extensions of 
ID schemes; (iii) triangular and tetrahedral meshes be used in 2D and 3D cases, respectively, so that the 
method is compatible to the simplest unstructured meshes and thus can be used to solve problems with 
complex geometries; and (iv) logical structures and approximation techniques used be as simple as possible, 
and special techniques that has only limited applicability and may cause undesirable side effects be avoided. 
Fortunately for the CESE development, as it turns out, the realization of the above lesser ideas (i)-(iv) 
follows effortlessly from that of the first two core ideas. 

The first model equation considered in the CESE development is the simple advection equation 


du du 

~di + a &c = 


( 1 . 1 ) 


where the advection speed a yf 0 is a constant. Let x\ = x, and X2 = t be considered as the coordinates 
of a two-dimensional Euclidean space E 2 . Then, because Eq. (1.1) can be expressed as V • h = 0 with 

h = f (au,u). Gauss’ divergence theorem in the space-time E 2 implies that Eq. (1.1) is the differential form 
of the integral conservation law 

( 1 l h ■ ds = 0 (1.2) 

Js(V ) 

As depicted in Fig. 1, here (i) S(V) is the boundary of an arbitrary space-time region V in E 2 , and (ii) 
ds = da n with da and n, respectively, being the area and the unit outward normal vector of a surface element 
on 5(V). Note that: (i) because h ■ ds is the space-time flux of h leaving the region V through the surface 
element ds, Eq. (1.2) simply states that the total space-time flux of h leaving V through S(V) vanishes; 
(ii) in E 2 , da is the length of a line segment on the simple closed curve S'(U); and (iii) all mathematical 
operations can be carried out as though E 2 were an ordinary two-dimensional Euclidean space. 

It is well known that a solution to Eq. (1.1) represents non-dissipative data propagation along its 
characteristic lines defined by dx/dt = a. Moreover, Eq. (1.1) is invariant under space-time inversion, i.e., it 
transforms back to itself if x and t are replaced by — x and —t, respectively. (In physics, space-time inversion 
invariance generally is referred to as PT invariance where P denotes a parity, i.e., mirror-image or spatial- 
reflection, operation while T denotes a time-reversal operation.) Thus a solution to Eq. (1.1) possesses the 
following properties: (i) it is completely determined by the data specified at an initial time level; (ii) its 
value at a space-time point has a finite domain of dependence (a point) at the initial time level; and (iii) the 
space-time inversion image of a solution to Eq. (1.1) is also a solution and vice versa. As such, in the initial 
CESE development, the focus is on the construction of an ideal core solver of Eq. (1.1) that enforces the 
conservation law Eq. (1.2) and also possesses properties similar to those of Eq. (1.1), i.e., it is a two-level, 
explicit, non-dissipative, and PT invariant solver. An in-depth account of this development and the resulting 
“a” scheme is given in [71]. As it turns out, the 2nd-orcler accurate a scheme (i) has a space-time stencil 
formed by one mesh point at the upper time level and two mesh points at the lower time level; and (ii) it 
is neutrally stable if v 2 < 1 where v = ciAt./Ax. Also, at each space-time mesh point ( j,n ), the a scheme 
is associated with two independent mesh variables uj and (uj)" (the numerical analogues of u and du/dx, 
respectively) and two equations. 

Until recently, with one exception (a three-level and 3rd-order accurate scheme reported on p. 80 of [1]), 
all CESE solvers of Eq. (1.1) are two-level and 2nd-order accurate extensions of the a scheme. To initiate 
a systematic CESE development of high-order schemes, in this paper we describe a new 4th-order accurate, 
space-time flux conserving, and neutrally stable CESE solver of Eq. (1.1). As will be shown, the space-time 
stencil of this two-level explicit scheme is formed by one point at the upper time level and three points at 
the lower time level. Because it is associated with three independent mesh variables u”, (« x )" and (u xx )™ 
(the numerical analogues of u , du/dx, and d 2 u/dx 2 , respectively) and three equations at each mesh point, 
hereafter the new scheme is referred to as the a(3) scheme. 
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The rest of the paper is organized as follows: In sec. 2, it is explained how the concepts of flux conserva- 
tion and PT invariance along with a requirement to minimize the truncation error of its principal component 
equation uniquely define the a(3) scheme. Also, for the a(3) scheme, we present (i) several of its equivalent 
forms; (ii) the truncation errors of its three component equations; and (iii) a set of PT-invariance induced 
algebraic relations involving the coefficient matrices of its component equations. 

A von Neumann analysis for the a(3) scheme is presented in Sec. 3. Specifically, we provide rather 
rigorous and thorough discussions on the properties of the 3x3 amplification matrix and its eigenvalues 
(i.e., the amplification factors). In particular, it is proved that: (i) the a(3) scheme must be neutrally stable 
if it is stable; (ii) all three amplification factors are of unit magnitude for all phase angles 9 of the Fourier 
modes considered in the von Neumann analysis if and only if \v\ < 1/2 ( v = aAt/Ax ); (iii) the a(3) scheme 
is stable if and only if \v\ < 1/2; and (iv) the a(3) scheme is linearly unstable (in a sense to be defined) if 

M = 1/2- 

In addition to numerically verifying the theoretical predictions made in Sec. 3, in Sec. 4 it is shown that 
the a(3) scheme generally is (i) 4th-order accurate for the mesh variables u” and (u x )"; and (ii) 2nd-order 
accurate for {u xx )^ . However, as predicted from theoretical considerations, in some exceptional cases the 
scheme can achieve perfect accuracy aside from round-off errors. Moreover, it is shown that the phase errors 
of the principal amplification factor of the a(3) scheme are 0{9 A ) if \v\ < 1/2, a sharp reduction from those 
of the dual a scheme [71] which are 0{6 2 ) if \v\ < 1. 

Conclusions and discussions are given in Sec. 5. Finally, several theorems and trigonometric identities 
used in Secs. 2 and 3 are proved in Appendices A and B while the three Fortran codes from which the 
numerical results presented in Sec. 4 are generated are listed in Appendices C, D, and E. 
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2. The a(3) scheme 

To proceed, consider the set fi of space-time mesh points (j, n) (marked by dots and crosses in Fig. 2(a)) 


where 

ft d = {(j,n)\j,n = 

0. ±1, ±2, ±3, . . .} 

(2.1) 

We have 

H = H 

1 U H 2 

(2.2) 

where Hi and H 2 are two disjoint sets defined by 



fti= f {(j,n)|j,n = 0,±l,±2,±3,. 

. . , and (j + n ) is an odd integer} 

(2.3) 

ft 2 = f {(j, n)\j, n = 0, ±1, ±2, ±3, . 

. . , and (j + n) is an even integer} 

(2.4) 


In Fig. 2(a), the mesh points G Hi are marked by dots while those £ fi 2 are marked by crosses. Hereafter 
O 2 is referred to as the complement set of Hi and vice versa. Obviously each of Hi and H 2 represents a set 
of space-time staggered mesh points. 

Each (j, n) G H is associated with (i) a solution element (SE), denoted by SE(j, n) (see Fig. 2(b) where 
(j,n) £ Hi is assumed), and (ii) two conservation elements (CEs), denoted by CE_(j, n) and CE+(j, n) (see 
Figs. 2(c) and 2(d) where (j,n) € Hi is assumed), respectively. Each SE is the interior of a space-time 
region that includes a horizontal line segment, a vertical line segment, and their immediate neighborhood. 
On the other hand, each CE is a rectangular space-time region. Hereafter, (i) SEs or CEs associated with 
mesh points £ Hi (G H 2 ) may be referred to simply as SEs or CEs associated with Hi (H 2 ). 

As a preliminary for the following development, note that (see Figs. 2(a)-(d)): 

(a) Two CEs which are associated with two mesh points, one G Hi while another G H 2 may occupy the 
same space-time region. As an example, (i) CE_(j, n) and CE + (_j — 1 ,n) occupy the same space-time 
region; and (ii) (j,n) £ Hi <*=> (j — 1, n) G H 2 . Hereafter the symbol “<t=>” is used as a shorthand for the 
statement “if and only if” . 

(b) A pair of diagonally opposite vertices of a CE both belong to the same set Hi or H 2 while another pair 
both belong to the complement set. As an example, points A and C belong to Hi while points B and 
D belong to fi 2 . 

(c) The CEs associated with each of Hi and H 2 by themselves are nonoverlapping and can fill the space-time 
E 2 . 

(d) Among the line segments forming the boundary of the same space-time region occupied by both 
CE_(j, n) and CE+(j — l,n), (i) AB and AD C SE (j, n); (ii) CB and CD C SE(j — l,n— 1); (iii) BA 
and BC C SE(j — l,n); and (iv) DA and DC C SE (j,n— 1). Because AB and BA represent the same 
line segment, one can see that any line segment on this boundary is a subset of two SEs with one of 
them being associated with Hi and another associated with H 2 . Hereafter, this ambiguity is removed 
by the following SE designation rule: any line segment designated as a boundary of a CE associated 
with Hi (H 2 ) is designated as a subset of a SE associated with Hi (H 2 ). As an example, if AB, AD, 
CB, and CD are designated as boundaries of CE_(j, n), then because points A and C belong to Hi, 
the above rule implies that: (i) both AB and AD are designated as subsets of SE(j, n); and (ii) both 
CB and CD are designated as subsets of SE(j — 1, n — 1). On the other hand, if BA, BC, DA, and DC 
are designated as boundaries of CE + (j — l,n), then: (i) both BA and BC are designated as subsets of 
SE(j — 1 ,n); and (ii) both DA and DC are designated as subsets of SE(j, n — 1). 

Let ( x,t ) G SE(j, n). Then Eqs. (1.1) and (1.2) will be simulated numerically assuming that u(x,t) and 
h(x,t), respectively, are approximated by 

u*(x,t . ; j, n) = f u" + ( u x )j {x-x j ) + {u t ) 1 -{t-t n )+^ (u xx )™ (x-Xj) 2 + (u xt )" (x-Xj)(t-t n ) + ^ ( u tt )j (t-t n ) 2 

(2.5) 
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and 


h*(x,t\j,n ) d = (au*(x,t-,j,n), u*(x,t;j,n )) (2.6) 

Note that: (i) w", (u x )", (ut)", (u xx ) r -, (w x t)" , and (utt)j are constants in SE(j, n), and the numerical 
analogues of the values of u, du/dx, du/dt , d 2 u/dx 2 , d 2 u/dxdt 1 and d 2 u/dt 2 at the mesh point ( j,n ), 
respectively; (ii) (xj,t n ) are the coordinates of the mesh point (j, n) where = jAa; and f” = nAt\ (iii) 
u*(x,t;j,ri) represents a 2nd-order Taylor’s approximation of u: and (iv) Eq. (2.6) is the numerical analogy 
of the definition ft = ( au,u ). 

For any (j, n) G O, let u = u*(x,t;j,n ) satisfy Eq. (1.1) for all (x,t) G SE (j,n). Then one has 

{ut)j = -a(u x )j, {u xt }] = -a(u xx )j , and (u tt )] = a 2 (« xx )"i (j,n)ef2 (2.7) 

Substituting Eq. (2.7) into Eq. (2.5), one has 

u*(x,t\j,n ) = u] + (u x )j [(a; - Xj) -a(t- t n )\ + ^( m xx )” [(a: - Xj) -a(t- t n )]~, (, j,n ) G Q (2.8) 


i.e., u", (uj,)", and ( u xx )" are the only independent mesh variables associated with (j, n). 

With the above preliminaries, next we derive the flux conservation relations that underline the a(3) 
scheme. 

2.1. Flux conservation relations 

Let the flux of h* conserve over all CEs, i.e., 


and 


h* ■ ds = 0, 

0») 

(j, n) G O 

(2.9) 

ft* ■ ds = 0, 

(j, n) G O 

(2.10) 


S(C£+0») 


Because (i) with respect to CE_(j, n), the outward unit normal vectors n at AB , AH, CD, and CB are 
(0, 1), (1,0), (0, —1), and (—1,0), respectively; and (ii) with respect to CE + (j, n), the vectors n at AF, AD , 
ED , and EF are (0,1), (—1,0), (0,-1), and (1,0), respectively, by using (i) the definitions given following 
Eq. (1.2), (ii) the above SE designation rule, and (iii) Eqs. (2.6) and (2.8), it can be shown that Eqs. (2.9) 
and (2.10) are equivalent to 


(1 + *) 
and 
(1-*) 


, x 2(1 — v + v 2 ) 
u — (1 — v)u x H u xx 


= {l + v) 


J j 


\ 2(1 + v + 

U + (1 + V)u x H u x 


respectively. Here: (i) v = f aAt/Ax is the Courant number; (ii) 

def AX 


= (1 -v) 


J 3 


x 2(1 — v + v 2 ) 

U + (1 - v)u x H u xx 


. x 2(1 + v + v 2 ) 
u - (1 + v)u x H u x 


(«*)"= and («ss)r= ^-r-( u ~)5 


def (Ax) 2 

4~ 


- n— 1 

? 

(j, n) G n 

Jj-1 

(2.11) 

- n— 1 

5 

(j, «) G 

■J j+i 

(2.12) 


(2.13) 


and (iii) to simplify notation, in the above and hereafter we adopt a convention that can be explained using 
an expression on the left side of Eq. (2.12) as an example, i.e., 


x 2(1 + v + is 2 ) 

U + (1 + v )u x H Uxx 


— + (1 + J/ ')( u x)? + 


2(1 + is + is 2 ) 


(Uxx) 1 - 


J J 
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(a) 


At this juncture, note that: 
Because 

du ax du 

dx 2 dx 


and 


d 2 u 
dx 2 


(a.t) 2 d 2 u 
4 dx 2, 


if 


_ def % 

X = a^/2 


the normalized parameters (u^)" and (uxx)”> respectively, can be interpreted as the numerical analogues 
of the values at (j, n ) of the first and second derivatives of u with respect to the normalized coordinate 

x. 

(b) By definition, points B and D depicted in Fig. 2(c) do not belong to either SE(j, n) or SE(j — 1, n — 1). 
This fact, however, does not pose a problem for flux evaluation over S(CE_(j,ri)) because the values 
of h* at isolated points do not contribute to the flux of h* over a finite line segment. Similarly, the fact 
that points D and F depicted in Fig. 2(d) do not belong to SE(j, n) and SE(j + 1, n — 1) does not pose 
a problem for flux evaluation over S(CE + (j,n)). 

(c) According to the SE designation rule, each line segment such as AB depicted in Fig. 2(c) can be 
assigned with two different fluxes of h* , one is associated with Q -j (hereafter referred to as the fii-flux) 
and another associated with Cl 2 (hereafter referred to as the f^-fhix). As such, among those local 
conservation relations Eqs. (2.9) and (2.10), those associated with (j,ri) £ fli are completely decoupled 
from those associated with (j, n) € SI 2 - Because Eqs. (2.9) and (2.10) are equivalent to Eqs. (2.11) and 
(2.12), respectively, it follows that each of the two systems of equations defined by Eqs. (2.11) and (2.12) 
is formed by two decoupled subsystems, one is associated with Q 1 while another associated with fl- 2 . 

(d) Moreover, because (i) the vector h* at any interface separating two neighboring CEs associated with the 

same set 12 1 (O 2 ) is evaluated using the information from the same SE, and (ii) the unit outward normal 
vector on the surface element pointing outward from one of these two neighboring CEs is exactly the 
negative of that pointing outward from another CE, one concludes that the flux leaving one of these CEs 
through the interface is the negative of that leaving another CE through the same interface. Due to this 
interface flux cancelation and the fact that the CEs associated with each of fii and 0,2 by themselves 
are nonoverlapping and can fill the space-time E 2 , the local conservation relations Eqs. (2.9) and (2.10) 
associated with (j,ri) £ (( j,n ) £ fl 2 ) lead to a global conservation relation, i.e., the total Q-] - (f l 2 -) 

flux of h* leaving the boundary of any space-time region that is the union of any combination of CEs 
associated with the same set 12 1 (Cl 2 ) vanishes. 


Let 1 — v 2 0, i.e. 1 + v ^ 0 and 1 — v ^ 0. Then Eqs. (2.11) and (2.12) reduce to 


, \ 2(1 — v + v 2 ) 

u - (1 - r)«i -I 


, x 2(1 — v + v 2 ) 

U + (1 - u)Ux H Uxx 


, (j,n)efi (2.14) 


J j-i 


and 


, \ 2(1 + v + v 2 ) 1" 

u + (1 + ^)n® H ^ Uxx 


n x 2(1 + v + v 2 ) - 1 n ~ 1 

u - (1 + U)Ux -\ 5 


, (j,n) £ ft (2.15) 


3 + 1 


respectively. Obviously, each of the two systems of equations defined by Eqs. (2.14) and (2.15) is also formed 
by two decoupled subsystems. Moreover, each component equation in Eq. (2.14) represents a stronger 
condition than the corresponding equation in Eq. (2.11) in the sense that the former implies the latter 
for any given v while the latter implies the former only if an extra condition (i.e., v ^ —1 for this case) 
is imposed. Similarly, each component equation in Eq. (2.15) represents a stronger condition than the 
corresponding equation in Eq. (2.12). These stronger conditions will be used in the construction of the a( 3) 
scheme. 

As a preliminary to a later development, next we will take a side tour and introduce the concept of 
invariance under space-time inversion. 

2.2. Invariance under space-time inversion 
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(2.16) 


Let u = u(x,t) be a solution to Eq. (1.1) in the domain — oo < x,t < +oo, i.e., 
du(x, t) du(x, t) 


dt 


dx 


= 0, 


— 00 < X, t < +00 


Let 

and 

Then (i) Eq. (2.16) +> 
and (ii) 

Thus Eq. (2.16) <+> 


x' = f — x and t! = f —t 


u(x , t) = f u(—x, —t) 


du{x',t') du(x',t') 


d d 

— — — — — — - and 
dt' dt 


dii(x, t) | du(x, t) _ q 


—oo < x ' , t! < +oo 


d_ _ _d_ 
dx' dx 

— 00 < X. t < +00 


(2.17) 

(2.18) 

(2.19) 

(2.20) 
(2.21) 


dt dx 

In other words, if u = u(x,t) is a solution to Eq. (1.1), so must be u = u(x,t ) and vice versa. Because the 
one-to-one mapping 

(x, t) <-> (— x, — f), —oo<x,t<+oo (2.22) 

represents a space-time inversion ( PT ) operation, hereafter (i) a pair of functions such as u and u will be 
referred to as the PT images of each other; and (ii) a partial differential equation (PDE) such as Eq. (1.1) 
is said to be PT invariant if the PT image of a solution is also a solution and vice versa. 

Next let 

u (k ’ e \x,t ) d = 9 d J^ and u {k ’ e \x,t) d = ® > -oo < x, t < +oo; k, £ = 0, 1, 2, . . . 

(2.23) 

Then, with the aid of the chain rule, Eqs. (2.17), (2.18), and (2.23) imply that 


)(M) 


(x, t ) = 


d k+e u(-x,-t) = k+e d k+e u(x',t') 


dx k dt l dx ,k dt fi — oo < x, £ < +oo; k,£ = 0,1,2,. 

= (-1 ) k+l u^\x\t') = (-1 ) fc+ V M) (-+-*) 


i.e., 


u^(x,t) = 


_ J u^ k ’^ (— x , — t ) if ( k + £) is even 


(2.24) 


(2.25) 


—u( k ’ e \— x, — t ) if ( k + £) is odd 

According to Eq. (2.23), u^ 0,0 ^ = u and til 0,0 ) = u. Thus Eq. (2.18) is a special case of Eq. (2.24) with 
k = £ = 0. 

In the following, the concept of PT invariance will be introduced for the a(3) scheme. As a preliminary, 
note that: (i) 

{j, n) <-> (-j, -n) (2.26) 

is the numerical analogue of the PT mapping Eq. (2.22); and (ii) (u x )™, ( u t )", {u xx )j, («xt)"> and (u tt )" 
are the numerical analogues of the values of u, du/dx, du/dt , d 2 u/dx 2 , d 2 u/dxdt 1 and d 2 u/dt 2 , at the mesh 
point ( j,n ), respectively. Thus, motivated by Eq. (2.25), the one-to-one mapping 


< - «=?; («,)• - -(«*)!?; K)" - -(«*)= 


3 —J 

>j‘ ‘ 


/j 

-j ) ^ XlJ—jl \^ZZJj 


— n 
3 




(2.27) 
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is taken as the numerical analogue of the one-to-one mapping 


’(x,t) <-> —oo<x,t<+oo-, k,£ = 0,1, 2, 3 

For the independent mesh variables, by using Eq. (2.13), Eq. (2.27) reduces to 


-(«*)=? 

(«**):" 


Eq. (2.29) can be expressed as 


O', n) € O 


q{j,n) <-> U q(—j, —n), 


0, «) G 


where 


9lj>)= (4? 
V («**)? 


0, n) G fi 


/ 1 0 0 \ 

U d = 0 -1 0 (2.32) 

Vo 0 \) 

The matrix U is unitary. In fact it is a real matrix with 

U = U - 1 (2.33) 

Hereafter (i) M -1 denotes the inverse of any nonsingular square matrix M; (ii) for each ( j,n ), Uq(—j, —n) 
is referred to as the PT image of q(j, n); and (iii) the set formed by Uq(—j, —n), ( j , n) G H is also referred to 
as the image of the set formed by q(j, n), ( j , n) G O. According to Eq. (2.33), q(j, n ) = UUq(—(—j), — (— n)). 
Thus q{j,n) is the PT image of Uq(—j,—n ) as an individual (j, n) or as the set defined over fi. In the 
following, we will show that by itself each of the four subsystems of equations associated with Eqs. (2.14) 
and (2.15) is PT invariant, i.e. , the subsystem maps onto an equivalent subsystem under the mapping 
Eq. (2.29). 

As an example, consider the subsystem of equations formed by the component equations associated with 
Hi in Eq. (2.14). Let it be denoted as Eq. (2.14a). Under the mapping Eq. (2.29), Eq. (2.14a) maps onto 


, . 2(1 — v + v 2 ) 1 " [ 2(1 — v + v 2 ) 1 ^ 

u + (1 - V)Ux -\ Uxx = U - (1 - V)Ux J \ Uxx 

6 1-3 L 6 J -0-1) 


, (j,n)GOi (2.34) 


At this juncture, note that, in addition to changing the sign of each Ux, mapping Eq. (2.29) requires that the 
upper and lower indices j, n, j — 1, and n— 1 in Eq. (2.14a) be replaced by their negatives, respectively. This 
is different from simply replacing the symbols j and n everywhere with —j and —n, respectively . Moreover, 
to simplify argument, hereafter system B is referred to as the PT image of system AHA maps onto B under 
the mapping Eq. (2.29), e.g., the subsystem Eq. (2.34) is the PT image of Eq. (2.14a). Let 

j* == 1 — j and n* d = 1 — n, ( j , n) G Hi (2.35) 

Then, by using the fact that ( j * + n*) + (j + n) = 2 and therefore ( j*,n *) G fii <t=> (j, n) G Hi, Eq. (2.34) 
can be cast into the form 


^ i 2(1 — v + v 2 ) 2(1 — v + v 2 ) 1" 

U — (1 — v)Ux -I Uxx = U + (1 - v)Ux H Uxx 

6 Jj* L d Jj*-1 


, (i*,n*) G Hi (2.36) 
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By comparing Eqs. (2.14a) and (2.36), one can see that the subsystem Eq. (2.14a) is identical to its PT 
image Eq. (2.34) (which is identical to Eq. (2.36)). Thus, under the mapping Eq. (2.29), Eq. (2.14a) maps 
onto itself, i.e. , the subsystem Eq. (2.14a) is PT invariant. QED. 

The PT invariance of another three subsystems associated with Eqs. (2.14) and (2.15) can be established 
in a similar manner. As such the system formed by all component equations in each of Eqs. (2.14) and (2.15) 
is PT invariant. 

The three mesh variables at any (j,n) £ £2 are linked to those at (j — l,n — 1) and (j + l,n — 1) by 
two component equations in Eqs. (2.14) and (2.15), respectively. In order that the three mesh variables at 
(j, n) can be determined in terms of those mesh variables at the (n — l)th time level, in the next subsection 
we introduce an extra PT invariant condition that links the mesh variables at (j, n) with those at the mesh 
point (j, n — 1). 

2.3. A family of PT invariant solvers 

Consider the following system of equations: 


[u + au s + fiuxx}™ = [u - au s + /3uxx}™ \ (j,n)e£2 (2.37) 

where a and (3 are parameters independent of (j, n). By definition, (j,n) £ fii (f^) +> ( j,n — 1) £ O 2 (£2i). 
Thus, unlike Eqs. (2.14) and (2.15), the mesh variables associated with £2i are linked to those associated with 
£2 2 through Eq. (2.37). However, as will be shown, like a subsystem associated with Eq. (2.14) or Eq. (2.15), 
the system of equations Eq. (2.37) is PT invariant for any pair of a and /3. 

The PT image of the system Eq. (2.37) is 


[u - aux + f3u m \_ " = [u + au s + Pux^J™ 1} , ( j , n) £ £2 (2.38) 

Let 

j' d = —j and n = f 1 — n, ( j,n)£fl (2.39) 

Then because ( j',n ') £ £2 4$ ( j,n ) £ f 2, Eq. (2.38) can be cast into the form 

[u + aUx + flUxx]y =[u- aux + fe]", -1 , (j', n') £ £2 (2.40) 

By comparing Eqs. (2.37) and (2.40), one can see that the system Eq. (2.37) is identical to its PT image 
Eq. (2.38) (which is identical to Eq. (2.40)). Thus, under the mapping Eq. (2.29), Eq. (2.37) maps onto 
itself, i.e., the system Eq. (2.37) is PT invariant. QED. 

Because each of Eqs. (2.14) and (2.15) is PT invariant, one can see that, for any pair of a and /3, the 
system formed by Eqs. (2.14), (2.15), and (2.37) is PT invariant. 

Next, the three mesh variables at any (j, n) £ £2 will be solved in terms of those at (j — 1, n — 1), (j,n— 1) 
and ( j + l,n — 1) using Eqs. (2.14), (2.15), and (2.37). Let 

A d = ^(\ + av) -2/3 (2.41) 

and assume a/ 0. Then it can be shown that Eqs. (2.14), (2.15), and (2.37) <t=> 



u — (1 + v)ux 
u + (1 — v)ux 


2(1 + v + v 2 ) 
3 

2(1 — v + v 2 ) 
3 


- n— 1 

- j+1 

- n— 1 

- 3 - 1 


(j,n)£ 12 (2.42) 
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1 n— 1 


4 ^ 


(%)" = ^ [u - au x + Pu xx ] 


n— 1 

j 


1 r2(l - i/ + v 2 ) 


A L 
1 

A L 


~P 


, . 2(1 + ^ + i^ 2 ) 

M - (1 + H U xx 


Jj+1 

2(1 + i/ + i/ 2 ) ir , 2(1 -v + v 2 ) i"-i 

o J L o i j- 1 


(j,n)Gfi (2.43) 


and 


(«s $ )” = [w - au x + (3u xx 1 " 

J A ' 


n— 1 

i 


1 — i' + a 


A 

1 + z/ — a 


2(1 + v + is 2 ) 1 n 1 

u - (1 + i/)Ua H u xx 

3 J j + i 

2(1 — V + I/ 2 ) -in— 1 

M + (1 - I/)«® H U xx 

O ij - 1 


(j, n) G ft 


(2.44) 


For any pair of a and /3 with A / 0, Eqs. (2.42)-(2.44) represent a solver for Eq. (1.1). In the next 
subsection, we pick out the pair of a and /3 with which the solver will have the smallest truncation error 
(i.e., the highest order of truncation error) for Eq. (2.42). 

2.4. A study of truncation error 

Because, at each Eqs. (2.42)-(2.44) represent a system of three equations for three independent 

mesh variables, Eqs. (2.42)-(2.44) represent a numerical analogue of a system of three coupled partial differ- 
ential equations (PDEs) with three dependent variables. (Eq. (1.1) is one of these PDEs). As such, in the 
following study, three different symbols u, v , and w will be used to denote the analytical versions of it”, and 
the non-normalized variables ( u x )" and ( u xx )”, respectively. Specifically, let u(x,t), v(x,t), and w(x,t ) be 
functions having all the derivatives needed. Thus one can define 




dx 


and 


w(x, t) ^ w(x, t) 


d 2 u(x, t) 
dx 2 


Also, as an example, one can define 


/ gi+^u 
\dx e dt m 


def d e+ u . 

= ^ a n — (iAx,nA t.) 
dx e dt m u ’ ’ 


£, ?n = 0,1,2,... 


(2.45) 


(2.46) 


Next we will consider the “analytical” version of Eq. (2.42) which results from replacing (i) it”, ( u x )”, 
and (iijj)", respectively, with ft”, u”, and iZ>”, for each ( j,n ); and (ii) the index n with n + 1 everywhere. 
By using Eq. (2.13) and the fact that (j, n + 1) G ft <G> ( j , n) G ft, the analytical form can be expressed as 


/ \ n def 1 [ ~n+l 4 [~ Oi AX _ 


1 

' (2av 

rt \ / -i \ 

2a' 

A 

\~3~ ‘ 

- /?J (i - 1/) - 

T. 

1 

~ (2av 

\ 

2a' 

A 

\Y~ ‘ 

- /?] (1 + i/) + 

T. 

j,n) G ft; 

a/0 



P(ax ) 2 


{1 + v)ax„ 
u v - 


(1 — v)ax , 


(1 


i^ + ^ 2 )Aa’ 2 i n 
— w 


6 


J l+i 


(1 — + v 2 )ax 2 ] n 

c w 

6 Jj-i 


(2.47) 


= 0 


NASA/TM— 2008-215138 


11 



By applying Taylor’s formula, it can be shown that 


(ei)? d = 


/ du du\ 4(a — v)duAx / d 2 u 2 d 2 u\At 1 r2ai/ 3 


\dt 


dx 


3a dx At V dt 2 


dx 2 ) 2 


2v(v — a) d 2 u (ax) 2 (d 3 u 3 d 3 u^ (At) 2 a d 2 v (ax) 
3a dx 2 At 


a L 3 

3 


+ P(l-V 2 ) 


2a . 


I ~~(1 T V 2 -\~ U 4 ) — (3lS 3 

3a L 3 ’ J dx At 

| 4 „~, £&n.\ ( 1 r 2ctl/3 


V dt 3 <9x 3 / 6 3a dx 2 At 

dw (ax) 3 a( 1 + 4zA) — 3/3^ — 2i/ 3 d 3 ft (ax) 3 


/ d 4 u ^d 4 ii\ (At) 3 J_ 


Vdf 4 
1 


<9x 4 7 24 6 a L 3 


+ P(l-v 2 ) 


9a 

d 3 u (ax) 4 
(9x 3 At, 


dx 3 At 

(3 d 2 w (ax) 4 
6a dx 2 At, 


12a nr + (1+2i/2 ~ i/ 


2a , 




d 4 u (ax) 

dx 4 At 


4 ^ « 


0[(At) 4 


— (1 - V + u 2 ) + (3(1 - I/)J [0[(AX) 5 ]/At + 0 [(ax) 4 ] + 0[At(AX’) 3 ] 

^ ^-(1 + V + v 2 ) - (3(1 + v) 0[(AX) 5 ]/At + 0 [(ax) 4 ] + 0[At(Ax) 3 ] 


dv (ax) 2 
dx At 


(2.48) 


(j, n) € O ; A ± 0 


Note that (ei)" defined in Eq. (2.47) is normalized by the factor (1 /a t) so that the lowest-order terms in 
the above Taylor’s expansion contain the leading term (du/dt, + adu/dx) which is independent of At and 
ax. Also, in Eq. (2.48) a term is denoted by 0[(At) fl (ax)^ 2 ] if there exists a constant C > 0 and two fixed 
integers i\ > 0 and £2 > 0 such that the absolute value of this term < C(At.) (l (ax) < 2 for all sufficiently small 
At and ax. Note that, in determining the order of magnitude of a term such as O [(ax) 5 ] in Eq. (2.48), 
the parameters a and (3 are not assumed to be constants independent of At and Ax. In fact, to reduce the 
truncation error of the a(3) scheme, they will be chosen to be functions of v (see Eqs. (2.58)) and thus vary 
with the ratio At/ Ax. 

In the following, let u = u(x, t), v = v(x, t), and w = w(x, t) be a solution to the system of PDEs formed 
by Eq. (1.1) and 

du d 2 u 

v - ^ = 0 and w — jrw = 0 (2.49) 


i.e., 


du du 

m +a d^ = 0 ’ 


du 

u~ — = 0 , 
dx 


and 


d 2 u 

dx 2 


= 0 


(2.50) 


In other words, here the scheme Eqs. (2.42)-(2.44) is considered as a solver of the system of PDEs Eqs. (1.1) 
and (2.49). Eqs. (2.45) and (2.50) imply that 


Qi+mjj d e+m w 

dx e dt m 0 a dx e dt. m 

d 2 u 2 d 2 u _ / d d\ 
dt, 2 a dx 2 \dt a dx) 

d 3 u 3 d 3 u _ / d 2 d 2 2 
~di 3+a d// 3 = V dt, 2 ~ a dtd/r + ° 

d 4 u 2 d 2 \ f d 

dt, 4 a dx 4 \ dt 2 ° dx 2 ) \dt 


£,m = 0,1,2,... 


du du , \ 


d 2 \ f du du\ 
dx 2 ) \dt a dx ) 


d \ / du du 
~ a ~dx ) (dt +a £h 


= 0 



(2.51) 

(2.52) 

(2.53) 

(2.54) 
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Note that the first equation in Eq. (2.50), and Eqs. (2.52)-(2.54) are all special cases of 


flt+m [ f)k~ f) k qi' 

W + ^ ‘■‘“V = °' f.m = 0.1,2, ..<! *=1,2,3... 


With the hint provided by Eqs. (2.52)-(2.54), Eqs. (2.55) can be proved using the first equation in Eq. (2.50) 
and elementary algebra. 

By using Eqs. (2.46) and (2.51)-(2.54), one can see that (ei)" reduces to 


(ei)? = 


4 (a — v) du Ax 2 v{v — a) d 2 u ( ax ) 2 a(l + 4 v 2 ) — 2>f3v — 2i/ 3 d 3 u (ax) 3 

3a dx At 3a dx 2 At 9a dx 3 At 


\ 

3a 

1 

r2i/ 4 

12a 

l 3 




2av\i d 4 u (ax) 4 1 ” 


dx 4 At 


0[(At) 4 


— -^(1 — V + z/ 2 ) + (3(1 — u) 0 [(ax) 5 ~\/ At + 0 [(ax) 4 ] + 0[At(Ax) 3 ] 

t ^-(1 + V + v 2 ) - (3(1 + v) 0[(Ax) 5 ]/At + 0 [(ax) 4 ] + 0[At(Ax) 3 ] 

O', n)eil;A/0 


By definition, the expression on the right side of Eq. (2.56) represents the truncation error of Eqs. (2.42) 
if the scheme Eqs. (2.42)-(2.44) are considered as a solver of the system of PDEs Eqs. (1.1) and (2.49). 
Here the values of a and (3 will be chosen so that the truncation error will reach the highest order. From 
Eq. (2.56), one can see that the coefficients of the three lowest-order terms in the truncation error vanish if 

ot-v = 0 and a(l + Av 2 ) - 3/3v - 2v 3 = 0 (2.57) 

For the case v ^ 0, Eq. (2.57) +> 


a = v and (3 = 


1 + 2v 2 


Next the a(3) scheme will be defined as the special solver with a and (3 being chosen according to Eq. (2.58). 

2.5. The basic and forward marching forms of the a(3) scheme 
Assuming Eq. (2.58), Eqs. (2.37), (2.41)-(2.44) and (2.56) reduce to 


1 + 2u 2 ] n r 

U + VU X + “ U xx — U VU X 

6 J j L 


1 + 2v 2 


n O r 1 + 2v 2 in— 1 1 + 

Uj =2 [u- vux H Uxx \j 2 


A = 2/3 
— u — (1 + v)i 


2(1 + ^ + v 2 ) ] n - 1 

^ U'XX 

3 -lj+1 


u + (1 - v) 


. 2(1 — v + v 2 ) l"- 1 

v)u s + , U xx 

3 J j-i 


(j,n)eSl (2.61) 


( u s) 1 j — 2 12 [u — viix H <x xx \ j i 2 

1 + 2v \ 

2 — [«+(!- v ) 


1 + 2v 2 n n—1 1 — 2.V 

^ Uxx „• “f ^ 


1 + 2 v 
2 


, 2(1 — v + v 2 ) l" -1 

VjUx + r Uxx 

o -1 .7—1 


U—( 1 + v) 

- n—1 

\x 

J .7 —1 


2(1 + v + v 2 ) l"- 1 

U x + z Uxx 

3 j .7+1 
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/ \n q T | 1 T 2zz i n 1 

y^xx) j — 3 (a VUx T ~ Uxx\ j 


3 

2 L 


, \ ^(l + + Z^) 

U - (1 + u)Ux H Ji« 


- n— 1 

- i+i 


3 

2 L 


2(1 — + z^ 2 ) i n - 1 

u + (1 - v)Ux H ZtSx 

o Jj-1 


(j,n)£ n (2.63) 


and 


1 / d 4 u\ n |"(ax) 


(ei)j " 24 \ 9a: 4 J 


At 


+ 2a 2 At.(Ax ) 2 + a 4 (At ) 3 +0[(a<) 4 ] + 0[(ax) 4 ] 


+ 0[At(Ax) 3 ] + 0[(At) 2 (Ax)^ 


0[(ax) e 


(j,n)£Cl (2.64) 


At 


Note that: (i) the forms of the last four terms in Eq. (2.64) have been simplified using the definition 
v = aAt/ AX', and (ii) the expression on the right side of Eq. (2.64) represents the truncation error of 
Eq. (2.61) if the scheme formed by Eqs. (2.61)-(2.63) is considered as a solver of the system of PDEs 
Eqs. (1.1) and (2.49). 

Next we convert Eq. (2.62) into its analytical form by replacing (i) w™, ( u x )”, and (u xx )(j , respectively, 
with m™, v ", and w™, for each (j, n); and (ii) the index n with n+ 1 everywhere. By using (i) Eq. (2.13), (ii) 
v = a At I Ax, and (iii) the fact that ( j , n + 1) £ Cl <t=> (j, n) £ Cl, then after a normalization by the factor 1/2, 
the analytical form can be expressed as 


n def 1 ~ n+ 1 2,aAt 

j ~ 2 j (ax) 2 


(e 2 )" = -v; 1 * - 

— ( 

!ax V 

— ( 

!ax V 


^ 2aAt\ 

2ax V ax ) L 

„ 2aAt\ 

, 1 H 

2ax V ax / L 


aAt _ (ax ) 2 + 2a 2 (At) 2 i 71 

n -yr v + i2 % 

Ax + aAt_ (ax ) 2 + aAt.Ax + a 2 (At ) 2 i n 

u v + - —w 

2 6 J j + 1 

ax — aAt _ (ax ) 2 — aAt.Ax + a 2 (At ) 2 i n 

u H v + - — — — w = 0 

2 6 Jj-i 


(j, n) £ Cl (2.65) 


Similarly, the analytical form of Eq. (2.63) can be expressed as 

def 1 -,_li 2 

3 

1 


(e 3 )? = g< +1 + 


(ax ) 2 


aAt _ (ax) 2 + 2a 2 (At) 2 _ 

“-T- " + 


(ax) 2 

1 

(ax) 2 


u — 


ax + aAt _ (ax ) 2 + aAt ax + a 2 (At ) 2 


6 


-w 


3 + 1 


Ax — aAt „ (ax) 2 — aAt.Ax + a 2 (At ) 2 _ 

u+ v+ ^ 

2 6 Jj-i 


= 0 


(j, n) £ Cl (2.66) 


By using Taylor’s formula and Eq. (2.45), Eqs. (2.65) and (2.66) imply that 


(e 2 ) 


n 

j 


f. 

r dv 

dv d , 

f du 

du\i 

< V + 

— 

— ~ — ( 

— + a — 

. dt 

dx dx ' 

< dt 

dxJ\ 


d 2 v (At ) 2 
dt? 4 


d 2 v [(ax ) 2 — 2a 2 (At) 2 
dx 2 4 


dw a 2 (At ) 2 - (ax ) 2 d td 2 u 2 d 2 u\(At ) 2 (l + v 2 )d 3 u, 2 

+ dx 6 dx\dt 2 a dx 2 ) 4 12 dx 3 X 

+ 0[(At) 3 ] +0[(At) 2 Ax] +0 [a1(ax) 2 ] +0[(ax) 3 ] (j, n) £ Cl 


(2.67) 
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and 


(dv 

ra 2 

fdu 

du\ 

1 r, \ 


.dx 2 

y dt 

+ a d^J 


dt 


dw 
1 dx 


d 2 v 
dx 2 


At 

~6 


ra 2 

/'d 2 u 

2 d 2 u^ 
— o 

d 2 w 
H 

2 d 2 w I 
— 2a 2 

(Af) 2 ( 

fd 3 V 

a 2 m 

(AX) 2 

Lax 2 ' 

\dt 2 

dx 2 ) 

1 dt 2 

dx 2 \ 

12 

lax 3 

dx 2 ) 

1 6 


(j,n)€fl (2.68) 


(ax) 2 > + 0[(Af) 3 ] + 0[(At) 2 Ax + 0[At(Ax) 2 ] + 0 [(ax) 3 


(l + v 2 ) d 4 it 
12 dx 4 

Assuming Eqs. (2.51) and (2.55), [ei )" and (ea)” reduce to 
( d 3 u\ n [(ax) 2 + a 2 (At) 2 ] 


V<9x 3 / ■ 


W = - ( ^ 

and 


12 


+ 0[(At) 3 ] +0 [(At) 2 Ax] + 0[At(Ax) 2 ] + 0[(ax) 3 ] (j,ri)£fl (2.69) 


= - (£)” 1(Al)2 ) ° W1 +0[(At)»] +0[(Atf Ax] + 0 [a<(a»)»] +0[(AX) 3 


(j» € fi (2.70) 


respectively. 

Hereafter, for any v, let the system of equations defined by Eqs. (2.14), (2.15), and (2.59) be referred to 
as the basic form of the a(3) scheme while that defined by Eqs. (2.61)-(2.63) be referred to as the forward 
marching form of the a(3) scheme. Because (i) Eqs. (2.14), (2.15), and (2.37) <t=> Eqs. (2.42)-(2.44) if A ^ 0, 
and (ii) Eqs. (2.59)-(2.63) are special cases of Eqs. (2.37), and (2.41)-(2.44), respectively, the basic form of 
the a(3) scheme its forward marching form. Thus the essential conditions represented by these or other 
equivalent forms may be referred to simply as the a (3) scheme. 

With the above definitions, the expressions on the right sides of Eqs. (2.64), (2.69) and (2.70) represent 
the truncation errors of Eqs. (2.61)-(2.63), respectively, if the forward marching form of the a(3) scheme is 
considered as a solver of the system of PDEs Eqs. (1.1) and (2.49). According to Eqs. (2.69) and (2.70), 
(e 2 )” — » 0 and (e^)™ — > 0 as At, Ax — > 0, regardless how At and ax are related when At, Ax — > 0. On 
the other hand, Eq. (2.64) implies that (ei)” — > 0 as At, ax — > 0 only if the mesh refinement procedure is 
subjected to the condition 


(ax ) 4 

At, 


0 as At, ax — > 0 


(2.71) 


Thus the a(3) scheme is consistent with the system of PDEs Eqs. (1.1) and (2.49) if and only if Eq. (2.71) 
is satisfied. 

At this juncture, we offer the following remarks: 

(a) Let At / ax be held as constant as At, ax — » 0. Then for this mesh refinement procedure, Eqs. (2.64), 
(2.69), and (2.70) imply that the truncation errors for Eqs. (2.61)-(2.63), respectively, are third order, 
second order, and second order in At and ax. On the other hand, according to the numerical results 
presented in Sec. 4, the a(3) scheme generally is 4th order in accuracy for both u” and (w x )" while only 
2nd order in accuracy for («„)”. Note that order of truncation error and order of accuracy represent 
total different concepts (see Secs. 5 and 6 in [1]). The former is a measure of how well an analytical 
solution satisfies the discrete scheme while the latter represents a measure of how well a solution to 
the discrete scheme approximates the corresponding analytical solution. Thus the numerical results 
presented in Sec. 4 do not contradict the conclusion reached here. 

(b) Because (i) each of the two decoupled subsystems in each of Eqs. (2.14) and (2.15) is PT invariant by 
itself, and (ii) the system Eq. (2.37) is also PT invariant if a and (3 are parameters independent of 
( j,n ), by the definition of PT invariance one can easily see that the basic form of the a(3) scheme is 
PT invariant. 
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(c) Let q(j, n) = q 0 (j, n), (j, n ) £ f!, be a solution to the basic form. Then, by substituting q(j, n) = q 0 (j , n) 
into the basic form, one obtains a system of identities involving <f Q (j, n), ( j,n ) £ fi. Due to the PT 
invariance of the basic form, the above system of identities is equivalent to that obtained by substituting 
q(j,n) = Uq 0 (—j,—n) into the basic form. As such q(j,n) = q 0 (j,n), ( j,n ) £ fl, represent a solution 
to the basic form <=> q(j,n) = Uq 0 (—j,~n), ( j,n ) £ fi, represent another solution to the basic form. In 
other words, the PT image of a solution to the basic form is also a solution and vice versa. Obviously 
this conclusion is valid for other PT invariant forms of the a(3) scheme. 

Next, the forward marching form Eqs. (2.61)-(2.62) will be cast into a matrix form. Let 


-> r \ def 

co (A) = 


-v J, c+(v) d = [ -(l + J') ), c_(V) d = 

(1 + 2i/ 2 )/3 / \ (2/3) (1 + v + v 2 ) 


1 

1-v 
(2/3) (1 — v + v 2 ) 


_ ... 2\ ^ f-(l + u)/2\ ^ ( —(1 — v)/2 \ 

d 0 (v) d = \2v , d+{v) d = (1 — 2v)/2 , d-{v) d = — (1 + 2v)/2 

-3/ V (3/2) J V (3/2) 


2 -2v (2/3) (1 + 2v 2 ) 

Qo(v) = f do{v) [c 0 (/)] ; = [ 2v —2v 2 (2/3)i/(l + 2v 2 ) 


-3 3^ — (1 + 2v 2 ) 

— (1 + v)/2 (1 + v) 2 /2 —{l + i/)(l + t/ + v 2 )/3 ' 

Q+(v) d = f d + {v) [c+iv)) 1 = [ (1 - 2u)/2 -(1 - 2i/)(l + v)/2 (1 - 2i/)(l + v + ^ 2 )/3 

3/2 — (3/2) (1 + v) l + v + v 2 


and 


— (1 — v)/2 — (1 — J/ ) 2 /2 -(l-^)(l-^ + ^ )/3 

Q-{v) = d-(v) [c_(i/)] = [ — (l + 2i/)/2 — (1 + 2i/)(l — i/)/2 -(1 + 2z/)(l — i/ + ^ 2 )/3 

3/2 (3/2) (1 — v) l-v + v 2 


(2.72) 

(2.73) 

(2.74) 

(2.75) 

(2.76) 


Hereafter c* denote the transpose of any column or row matrix c. By using Eqs. (2.31) and (2.74)-(2.76), 
the forward marching form can be cast into the matrix form: 


q(j,n) = Qo(v)q(j,n — 1) + Q + {v)q{j + l,n- 1) + Q_(v)q(j -l ,n- 1), (. j,n ) € Q (2.77) 

Here the reader is warned that the notations <5+(z/) and Q-(u) used in earlier CESE papers are now replaced 
by Q~(v) and Q+(v ), respectively. As such, the terms Q-(v)q(j+ 1, n— 1) and Q + (i/)(j — 1, n—T) in Eq. (3.48) 
of [71] appear here as Q+{v)q(j + 1 ,n— 1) and Q_(v)(j — l,n — 1), respectively. Also note that each of 
Qo( v )i Q+{ v )i anc l Q-W) is * n the form of dc * where c and d are 3x1 column matrix. Thus each is a matrix 
of rank one (see pp. 80-82 in [74]). Rank-one matrices are singular and have many interesting properties. 
As an example, the eigenvalues of Qo(v) are 0, 0, and [co(z')]* do(v) with do(v) being the eigenvector of the 
last eigenvalue. 

To facilitate the proof of the PT invariance of the forward marching form, first we will introduce some 
basic concept. Note that, for any set of variables xg, ye, t = 1, 2, the conditions 

xi + y\ = x 2 - ?/2 and x\ - yi = x 2 + 2/2 (2.78) 

x\ = X 2 and 2/1 = ~V 2 (2.79) 

Thus, the image of Eq. (2.78) under any one-to-one mapping 

[xt,ye) {x' t ,y'e), f = l,2 (2.80) 
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1 .( 3 ., 

x'l + Vi = x' 2 - y 2 and x\ - y[ = x 2 + y 2 (2.81) 

<t=> the image of Eq. (2.79) under the same mapping, i.e., 

x\ = x 2 and y[ = - y 2 (2.82) 

where the variables x 'f, and y' e , £ = 1,2, may or may not be related to xe,ye, £ = 1,2. Moreover, in case 
that these two sets of variables are related, the condition Eq. (2.78) (or its equivalent Eq. (2.79)) may or 
may not be equivalent to the condition Eq. (2.81) (or its equivalent Eq. (2.82)). If the mapping Eq. (2.80) 
is such that Eq. (2.78) the image under this mapping (i.e., Eq. (2.81)), then Eq. (2.79) (the equivalent of 
Eq. (2.78)) <t=> Eq. (2.82) (the equivalent of Eq. (2.81)). Eq. (2.80) with x' e = xe and y' e = ye, £ = 1,2, is an 
example of such mapping while Eq. (2.80) with x\ = ye and y' e = xe, £ = 1,2, is not. 

To prove the PT invariance of the forward marching form, Note that: (i) the basic form of the a(3) 
scheme <t=> its forward marching form for any choice of q(j , n), ( j , n) G Cl; and (ii) the PT images of the basic 
and forward marching forms, respectively, are obtained from the basic and forward marching forms through 
the mapping Eq. (2.30), i.e., through replacing q{j, n) in the basic form and the forward marching form with 
Uq(—j,—n), (j,n) G Cl. From the above observations and the illustration given in the last paragraph, one 
concludes that the PT image of the basic form yy that of the forward marching form. Because the basic 
form is PT invariant, i.e., the PT image of the basic form the basic form itself, Now we arrive at the 
conclusion that the forward marching form the basic form <=> the PT image of the basic form <t=> the PT 
image of the forward marching form. Thus the forward marching form -o- its PT image, i.e., the forward 
marching form is PT invariant. QED. 

With the above preliminaries, the backward marching form of the a(3) scheme will be developed in 
Sec. 2.6. 

2.6. The backward marching forms of the a(3) scheme 

The PT invariance of the forward marching form of the a(3) scheme implies that Eq. (2.77) <t=> its PT 
image, i.e., 


Uq{-j,-n) = Q 0 (u)Uq{-j,-n+l) + Q + (u)Uq(-j - 1, -n+ 1) + Q-{v)Uq(-j + l,-n + 1), (. j,n ) G Cl 

(2.83) 

Moreover, by multiplying Eq. (2.83) from left using the matrix U and using Eq. (2.33), one concludes that 
Eq. (2.83) «• 


q(-j, -n) = Qoiy)q(-j , -n+ 1) + Q-{v)q(-j - 1, -n+ 1) + Q+{v)q(-j + 1, -n+ 1), 


(j,n) G Cl (2.84) 


where 


Qoiy) = UQ 0 (u)U = 


( -(I + 1O/2 

Q-[v) d = UQ+(is)U= — (1 — 2v)/2 
\ 3/2 

and 


/— (1 — 10/2 
Q+{v) A =UQ-(v)U = (1 + 2i0/2 

V 3/2 


/ 2 2v (2/3) (1 + 2v 2 ) \ 

—2v -2v 2 — (2/3)i/(l + 2v 2 ) 

\ -3 -3 v — (1 + 2v 2 ) ) 

— (1 + v) 2 /2 -(l + u )(l + u + u 2)/3 

— (1 — 2^)(1 + v)/2 — (1 — 2^)(1 + v + v 2 )/3 

(3/2) (1 + v) 1 + v + v 2 


(\-vf!2 
— (1 + 2i/)(l — v)/2 
(3/2) (1 — u) 


— (1 — i/)(l — v + ^ 2 )/3 
(1 + 2u)(l — v + r ,2 )/3 
1 — v + v 2 


(2.85) 


(2.86) 


(2.87) 


By replacing the “dummy” indices — j and —n everywhere in Eq. (2.84) with j and n, respectively, one can 
see that the system Eq. (2.84) is identical to the system 

q{j, n) = Q 0 {v)q{j, n + 1) + Q + {v)q{j + 1, n + 1) + Q-(u)q(j - 1 ,n+ 1), ( j , n) G Cl (2.88) 
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Because the mesh variables at (j,n) can be determined in terms of those at (j — 1 , n + 1), ( j,n + 1), and 
(j + l,ra+l) using Eq. (2.88), hereafter Eq. (2.88) (which is equivalent to other forms of the a( 3) scheme) 
will be referred to as the backward marching form of the a(3) scheme. 

According to Eqs. (2.74) and (2.85), Qo(^) = Udo(u) [c'o(^)] i U. Because Udo(v) and [c'o(z')]* U are 3x1 
column matrix and 1x3 row matrix, respectively, Q o(y) is a rank-one matrix. Similarly, Q_(y) and Q + (u) 
are also rank-one matrices. 

Eq. (2.88) was derived using the PT invariance of the forward marching form of the a(3) scheme. 
Alternatively, it can also be derived from the basic form. To proceed, note that: (i) by replacing the indices 
j and n everywhere in Eq. (2.14) with j + 1 and n + 1 and using the fact that (j, n) G ft -4+ (j — 1, n — 1) G 
ft <+> (j + 1, n + 1) G ft, one can see that the system Eq. (2.14) is identical to the system 


, \ 2(1 — v + v 2 ) 

U + (1 - U)U S H U S x 


n ^ 2(1- v + v 2 ) in+1 

u ~ (1 - v)u s -\ Uxx 

0 + 1 


(j,n)efl (2.89) 


(ii) by replacing the indices j and n everywhere in Eq. (2.15) with j — 1 and n+1 and using the fact that 
(j. n) € ft <+> (j + 1, n — 1) G 0 <+> (j — 1, n + 1) G 0, one can see that the system Eq. (2.15) is identical to 
the system 

n+1 

, (j, n) G (2.90) 

j - 1 

and (iii) by replacing the index n everywhere in Eq. (2.59) with ?i+ 1 and using the fact that (j, n) G +> 
( j , n — 1) G O <+> ( j , n + 1) G fi, one can see that the system Eq. (2.59) is identical to the system 


\ 2(1 + v + v ) 

U - (1 + v)Ux H Ux 


/i \ 2(1 + v + v ) 

u + (1 + V)Ux ~\ Ux 


U — ISUx 


1 + 2tA 1 



U + VUx + 


1 + 2v 2 


-i n+1 



(j, n) G 


(2.91) 


As such the system Eqs. (2.89)-(2.91) are identical to Eqs. (2.14), (2.15), (2.59), respectively. 

For each ( j,n ) G ft, Eqs. (2.14), (2.15), and (2.59) form a linear system of three equations for the three 
mesh variables w", (%)", and («xs)"- Eqs. (2.89)-(2.91) form another system. Moreover, one can see that, 
under the mesh variable mapping 


q{j, n) <-> U q(j, n ) , q{j, n - 1) <-> q(j, n + 1) , 

q{j + 1,71 - 1) <-> Uq(j - l,n + 1), and q(j - 1, n - 1) q(j + 1, n + 1) 


Eqs. (2.89)-(2.91), respectively, are the images of Eqs. (2.14), (2.15), and (2.59) and vice versa. By using 
the concept introduced earlier in a discussion involving Eqs. (2.78)-(2.82), one concludes that the solution 
to Eqs. (2.89)-(2.91) must be the image of Eq. (2.77) (i.e., the solution to Eqs. (2.14), (2.15) and (2.59)) 
under the same mapping. In other words, the solution to Eqs. (2.89)-(2.91) is 


Uq(j,n) = Qo(v)Uq(j,n+ 1) + Q + {v)Uq{j - l,n + 1) + Q-{v)Uq(j + l,n + 1), (. j,n ) G ft (2.93) 

By multiplying Eq. (2.93) from left using the matrix U and using Eqs. (2.33) and (2.85)-(2.87), one has 
Eq. (2.88). QED. 

As a preliminary for the developments in Sec. 3, in the following, important algebraic relations involving 
Qo{v), Q+(is), Q-(v), Qo{v), Q+{v), and Q~{v) will be extracted from the PT invariance of the a(3) scheme. 

2.7. Algebraic relations associated with PT invariance 

Let (jo, n 0 ) G ft be any given fixed mesh point. Let q\j 0 , n 0 ), q(jo ± 1, n G ), and q(j 0 ± 2, n 0 ), respectively, 
be the arbitrary initial data specified at (j D , n Q ), (j 0 ± 1, n Q ), and (j Q ± 2, n 0 ), respectively. Let q(j 0 , no + 1), 
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and q(j Q ± l,n 0 + 1) be specified in terms of the mesh variables at the n Q th time level using the forward 
marching form Eq. (2.77), i.e., 


q(jo, n a + 1) = Qo{v) q{jo, n 0 ) + Q+{v) q{j a + 1, n a ) + Q-{v) q{jo ~ 1, n 0 ) (2.94) 

q(j 0 + 1, n 0 + 1) = QoM + 1, n 0 ) + Q+{v) q{jo + 2, n 0 ) + Q-(v) q(jo, n 0 ) (2.95) 

and 

q(jo - 1, n a + 1) = QoM - 1, n 0 ) + Q+{v) q{jo, n 0 ) + Q-(v) q(j Q - 2, n 0 ) (2.96) 

On the other hand, because Eq. (2.77) <t=> Eq. (2.88), q{j 0 ,n 0 + 1), q(j 0 ± 1 ,n 0 + 1), and q{j 0 ,n 0 ) must also 
be linked by Eq. (2.88), i.e., 

q(jo,n 0 ) = Q 0 (v) q(jo,n 0 + 1) + Q+(v) q(j Q + Mo + 1) + Q-(zy) q{j a -l,n 0 + 1) (2.97) 

Substituting Eqs. (2.94)-(2.96) into (2.97), one has 

[Qo(v)Qo(v) + Q+MQ-M + Q-MQ+M - i]q{jo,n 0 ) 

+ [Qo(v)Q+(v) + Q+(v)Qo(v)\q(jo + 1 ,n 0 ) + [Q 0 (v)Q-(u) + Q-WQoM^jo - Mo) (2-98) 
+ Q+MQ+M q(jo + 2, n 0 ) + Q-(y)Q_(y) q(j a - 2,n c ) = 0 


where I is the 3x3 identity matrix and 0 is the 3x1 null column matrix. 

Because Eq. (2.98) must be valid for any choice of q(j 0 ,n 0 ), q(j 0 ± l,n G ), and q(j 0 ± 2 ,n 0 ), the coeffi- 
cients matrices in front of these column matrices must be null identically, i.e., 

Qo(v)Qo(v) + Q+(v)Q~(v) + Q-(y)Q+(y) = / (2.99) 

Qo(v)Q+(v) + Q + (u)Q 0 (u) = 0 (2.100) 

Q 0 {v)Q-{v) + Q_{v)Q q {v) = 0 (2.101) 

Q+(v)Q+(v) = 0 (2.102) 

and 

Q_(i/)Q_(i/) = 0 (2.103) 

where 0 is the 3x3 null matrix. As an example, one can prove Eq. (2.99) by substituting into Eq. (2.98) 
each of the following sets of the initial data: (i) q(j 0 ± 1 ,n 0 ) = q(jo ± 2 ,n 0 ) = 0 and q[j 0 ,n 0 ) = (1,0,0)*, 
(ii) q{j 0 ± 1, n 0 ) = q{j Q ± 2, n Q ) = 0 and q(j a , n Q ) = (0, 1, 0)*, and (iii) q(j Q ± 1, n Q ) = q(jo ± 2, n Q ) = 0 and 
q(jo, n 0 ) = (0, 0, 1)*. 

Similarly, by substituting the backward marching relations 

q(jo, n a - 1) = Qo(v) q{jo, n 0 ) + Q+{v) q{jo + 1, n 0 ) + Q-(v) q(j 0 - 1, n 0 ) (2.104) 

q(jo + 1, - 1) = QoM q{jo + 1, rio) + Q-M) q(jo + 2, n Q ) + Q_(i/) <?(j 0 , n D ) (2.105) 

and 

Stio - Mo - 1) = Qo(v) q(j 0 - 1, n 0 ) + Q + (i/) q{j 0 ,n 0 ) + Q-(v) q(j a - 2,n c ) (2.106) 

into the forward marching relation 

q(jo, no) = Qo(y) q(j 0 , n Q - 1) + Q+(y) q(jo + 1, n Q - 1) + Q-(y) q{j a - 1, n a - 1) (2.107) 
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one has 


[Qo(v)Qo(v) + <2 + 0)<2_0) + Q-(u)Q + (v) - I]q(j 0 ,n 0 ) 

+ [Qo(v)Q+(v) + Q+{v)Qo{v)]q{j 0 + l,n 0 ) + [Qo(v)Q-(v) + Q-(v)Qo(v)]q{jo ~ l,n 0 ) (2.108) 

+ Q+(v)Q+(v) q{jo + 2 ,n 0 ) + Q-(v)Q-(y) q(j 0 - 2 ,n„) = 0 

Because Eq. (2.107) must be valid for any choice of q(j 0 , n 0 ), q(jo ± 1 , n G ), and q(j 0 ± 2 ,n 0 ), one concludes 
that 

QoMQoW + Q+^Q-M + Q_(i/)Q+H = I (2.109) 

QoMQ+M + Q+(v)Qo(v) = 0 (2.110) 

Qo(v)Q-(v) + Q-(v)Qo(v) = 0 (2.111) 

Q+(i/)Q+(i/) = 0 (2.112) 

and 

O-MO-M = o (2.H3) 

By using Eqs. (2.32) and (2.85)-(2.87), it can be shown that: (i) Eq. (2.99) <+> Eq. (2.109) <+> 

Qo(v)UQ 0 (v) + Q-(v)UQ-(y) + Q+{y)UQ+{y) = U (2.114) 

(ii) Eq. (2.100) +> Eq. (2.111) <+> 

Qo{y)UQ+{y) + Q-{y)UQ 0 {y) = 0 (2.115) 

(iii) Eq. (2.101) +> Eq. (2.110) +> 

Qo[y)UQ-[y) + Q+{v)UQ 0 {v) = 0 (2.116) 

(iv) Eq. (2.102) +> Eq. (2.113) +> 

Q-{y)UQ+{v)=Q (2.117) 

and (v) Eq. (2.103) +> Eq. (2.112) +> 

Q+{v)UQ-{v) = 0 (2.118) 

2.8. Other invariant properties and related algebraic relations 

By using Eqs. (2.32) and (2.74)-(2.76), one can show that 

Qo(-z') = UQ 0 (v)U, Q-(-is) = UQ + (u)U, and Q+(-v) = UQ-{v)U (2.119) 

By using Eqs. (2.85)-(2.87), one can also show that Eq. (2.119) <+> 

Q 0 {-v) = UQ 0 {v)U, Q-(-v) = UQ+(v)U, and Q+(-u) = UQ-{v)U (2.120) 


As will be shown, the above relations are linked with other invariant properties of the a(3) scheme. 

Let the advection speed a in Eq. (1.1) be considered as a variable parameter. Let u = u{x,t\a) be a 
solution to Eq. (1.1), in the domain — oo < x,t,a < +oo, i.e., 


Let 


du(x, t; a) ^ du(x, t;a ) _ ^ 


dt 


dx 


/ def ,/ def , , 

x = —x, t = t, and a = —a, 


—oo < x,t,a < +oo 


—oo < x,t,a < +oo 


(2.121) 

(2.122) 
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and 


(2.123) 


Then (i) Eq. (2.121) +> 


u(x, t; a) = f u(—x, t; 


-a) 


du( x', t'\ a') 
dt' 


,du{x',t';a') _ n 
“ dx' 


—oo < x' , t' , a ' < +oo 


and (ii) 

d d d d 

di’ = dt anc dp = -fa 

Thus one concludes that Eq. (2.121) <+> 


du{x, t; a) 
Ot 


„ dii(x,t;a ) _ n 
a — u, 

ox 


— 00 < X, t < +00 


(2.124) 


(2.125) 


(2.126) 


In other words, if u = u(x,t\a) is a solution to Eq. (1.1), so must be u = u{x,t\a ) and vice versa. Because 
the one-to-one mapping 


(x, t, a) <-> (— x, t, — a), ~oo < x,t,a < +oo (2.127) 

represents a combined spatial-reflection (parity) and advection direction reversal (ADR) operation, hereafter 
(i) a pair of functions such as u and u will be referred to as the PADR images of each other; and (ii) a PDE 
such as Eq. (1.1) is said to be PADR invariant if the PADR image of a solution is also a solution and vice 
versa. 

Because v = aAt/Ax, the numerical analogue of Eq. (2.127) is 

(j, n) (— j,n ) and v <-> —v (2.128) 

Motivated by an argument similar to that leads to Eq. (2.30) for PT mapping, the PADR mapping for the 
a (3) scheme is defined by 


q(j,ri) <-> Uq(—j,ri) and v — u, (2.129) 

Thus the PADR image Eq. (2.77) is 

Uq(-j,n) = Qo(—u)Uq(—j, n — 1) + Q + {-v)Uq{-j - l,n - 1) + Q-(-v)Uq(-j + 1 ,n- 1), (, j,n ) G ii 

(2.130) 

By using Eqs. (2.32) and (2.119), it can be shown that Eq. (2.130) <+> 

q(-j,n ) = Qo{v)q(-j,n - l) + Q-{v)q(-j - l,n - l) + Q+{v)q(-j + l,n - l), ( j,n ) G ii (2.131) 

By replacing the dummy index —j with j everywhere in Eq. (2.131) and using the fact that (— j,n ) G ii +> 
(j,n) G ii, one concludes that Eq. (2.131) <+> Eq. (2.77). Thus Eq. (2.77) is PADR invariant, i.e., it is 
equivalent to its PADR image. 

By exchanging the roles of x and t, one can define invariance under a combined time reversal and 
advection direction reversal operation. Because (i) this operation is equivalent to a PT operation followed 
by a PADR operation or vice versa, and (ii) Eq. (1.1) and the a(3) scheme are invariant under both PT 
and PADR operations, one concludes that Eq. (1.1) and the a(3) scheme are also invariant under the new 
operation. In fact, invariance of the a(3) scheme under this new operation can be proved using Eq. (2.120) 
(which is equivalent to Eq. (2.119)). 
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3. von Neumann analysis 


Let G(v, 9) be a 3 x 3 nonsingular complex matrix function of v and the phase angle 9 such that 

q(j, n) = e 1 ^ [G(u, 9)] n b , (j, n) € O; — oo < v < +oo; — 7r < 9 < x (i = y/— 1) (3-1) 

is a solution to Eq. (2.77) for all possible complex constant 3x1 column matrices b. Note that: (i) without 
any loss of generality, hereafter the domain of 9 is limited to —n < 9 < tt and, unless specified otherwise, 


def 


| [G(u, 9)} 1 1 for an integer n < 0, 


this domain will be assumed implicitly; and (ii) because [G{v, 

[G(v, 9)] n is not defined if n < 0 unless [G(v, 0)] 1 exists, i.e. , G(v, 9) is nonsingular. By substituting Eq. (3.1) 
into Eq. (2.77), one has 


[GM) - QoM - e ie Q+{y ) - e~ ie Q-{v )] [G(u, 9)] n b = 0, 


71 = 0, rhl, ±2, . . . 


(3.2) 


Because (i) [G(v, 0)]° = /, and (ii) b can be any complex constant 3x1 column matrix, Eq. (3.2) 


G(v, 9) = Q 0 (y) + e l6 Q+(v) + e~ M Q-{v) 


—iO/ 


(3.3) 


By definition, G(v,9) is the amplification matrix of the forward marching form of the a(3) scheme. Because 
Qo{v), Q+(y), and Q-(v) are real matrices, Eq. (3.3) implies that 


G{v,-9) = G{v,9) (3.4) 

Hereafter M denotes the complex conjugate of any matrix M. Also, with the aid of Eq. (2.119) and the 
relation U = U~ l , one has 

G{—v, 9) = UG(u, —9)U = UG{v, -9)U~ l (3.5) 

At this juncture, some comments on the dual roles played by the amplification matrix G{ v, 9) in de- 
termining the accuracy and the stability of the a(3) scheme are in order. Note that the von Neumann 
analysis represents essentially a rigorous discrete Fourier analysis performed for a Fourier mode of a solution 
to a linear marching scheme such as Eq. (2.77) assuming periodic spatial boundary conditions (see Sec. 4 
in [1]). The only difference between them is that the parameter 9 (which specifies a Fourier mode) in the 
von Neumann analysis can assume any value in the domain — 7 r < 9 < tt while that in the discrete Fourier 
analysis can only assume a set of K uniformly distributed discrete values within the domain — 7r < 9 < tt if 
the spatial domain is divided into K uniform mesh intervals. As such, the time evolution and therefore the 
accuracy of a Fourier mode of a solution to a linear scheme assuming periodic spatial boundary conditions 
can be determined using the corresponding amplification matrix (see Sec. 5 in [1]). Moreover, because a 
linear combination of solutions to a linear marching scheme is also a solution by itself, the time evolution 
of any Fourier mode of the round-off errors originally introduced during any marching step is also governed 
by the linear scheme and therefore it can also be determined using the amplification factor. As a result of 
this consideration and the facts that: (i) a scheme is stable if and only if these round-off errors will not 
be amplified without bound after many marching steps, and (ii) the spectrum of round-off errors generally 
covers all possible Fourier modes, i.e., all possible discrete values of 9, one concludes that, assuming periodic 
spatial boundary conditions, the a( 3) scheme is stable for a given v if and only if, for all possible K discrete 
values of 9 in the domain —n < 9 < tt, every element of the matrix [G(v, 9)] m remains bounded as the 
positive integer m — > +oo. Because the distribution of the allowed discrete values of 9 becomes very dense 
in the domain — 7r < 9 < tt for a large K , for simplicity, the K discrete values of 9 referred to in the above 
stability definition is replaced by all values of 9 in the domain — 7r < 9 < tt in Definition 1 of Sec. 3.9. 

In the following, we will show that the a(3) scheme must be neutrally stable when it is stable. 

3.1. Neutral stability of the a(3) scheme 
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By using Eqs. (2.114)-(2.118), one can show easily that 


U[Q 0 (v) + e ie Q-(is) + e~ i0 Q + {is)] U [Q 0 p) + e i0 Q+(is) + e~ i0 Q-(y)\ 

= [Qop) + e*Q+p) + e- i0 Q-{u)]U[Q o {v) + e <fl Q_p) + e~ i0 Q + (is)\ U = I 


(3.6) 


Thus Gp, 9) defined in Eq. (3.3) is nonsingular and its inverse is 

[G(i/,0)r* = U [Qo(v) + e l0 Q.(v) + e~ i0 Q + (y )] U (3.7) 

Indeed, with the aid of Eq. (2.85)-(2.87), Eq. (3.7) is what one obtains after substituting Eq. (3.1) into the 
backward marching form Eq. (2.88). Moreover, by using Eqs. (2.32), (3.3), (3.4), and (3.7), one has 

[Gp, f?)] -1 = UG{y, 0)[7 _1 (3.8) 


For each p, 9), let the three eigenvalues of Gp, 9) be denoted as oyp, 6), £ = 1,2, 3, respectively. They 
will be referred to as the amplification factors of the a(3) scheme. Because Gp, 9) is nonsingular, 


<7/p,0)^O, 1= 1,2,3 


(3.9) 


(see part (i) of Theorem 1 given below). Also, as will be shown, a ((is, 9), t = 1, 2, 3, satisfy the following set 
condition: 


crip, 6) ' cr 2 p, 6) 1 cr 3 p, 6) 


jcripP), cr 2 p, 6), a 3 p,0) | 


(3.10) 


Hereafter z denotes the complex conjugate of any complex number z. 

As a preliminary, first we introduce the following well-known matrix theorems: 

Theorem 1. Let A be a nonsingular N x N matrix with the eigenvalues Xe, £ = 1, 2, . . . , N. Then (i) 
Xi ^ 0, £ = 1, 2, . . . , N; and (ii) the eigenvalues of A^ 1 are 1/Xe, £ = 1, 2, . . . , N. 

Theorem 2. Let A be a N x N matrix with the eigenvalues Xg, £ = 1,2 , ,N. Then the eigenvalues 
of A , the complex conjugate of A, are Xg, £ = 1,2, ... ,N. 

Theorem 3. Let A and B be two similar N x N matrices, i.e., there exists a nonsingular TV x N matrix 
S so that B = S' -1 AS 1 . Then A and B have the same eigenvalues, counting multiplicity. 

Theorems 1 and 2 are proved in Appendix A, while Theorem 3 is proved on p. 45 of [76]. 


To prove Eq. (3.10), note that part (ii) of Theorem 1 implies that, for any (//, 6), the eigenvalues of 
[G(u, pp 1 are 1/<ji(v,0), £ = 1,2,3. Next, by using Theorems 2 and 3, and the fact that (t/ -1 ) -1 = U, 
one can see that the eigenvalues of the matrix on the right side of Eq. (3.8) are <Je(v,9), £ = 1,2,3. Thus 
Eq. (3.10) now is an immediate result of Eq. (3.8). QED. 

An immediate result of Eq. (3.10) is 


111 

<Ji(v,9) <t 2 PP) cr 3 pp) 


crip, 9) ■ cr 2 p, 9) • <t 3 (i/, 9) 


i.e., 

I <7i p, 0)| • p2p,0)l • p 3 p,0)| = 1 (3.11) 

As will be shown in Sec. 3.9, for any given is, a necessary condition for the stability of the a(3) scheme is 

|o^p, 0)| < 1, £ = 1,2,3 (3.12) 

Thus Eq. (3.11) implies that, for any given v, the a(3) scheme must be neutrally stable, i.e., 

|<7/p,0)| = l, £=1,2,3 (3.13) 
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if it is stable. As such, Eq. (3.8) does not imply neutral stability of the a( 3) scheme. However, it does imply 
that the scheme can only be neutrally stable (i.e., non-dissipative) if it is stable. Here we have reached this 
conclusion without using the explicit form of ay (y, 9), £ = 1,2, 3. 

At this juncture, note that one can obtain 

oy(-+ 9) = ae(v, -9) = oy(+ 8), 6= 1,2,3 (3-14) 

by using Eqs. (3.4) and (3.5) along with Theorems 2 and 3. 

Eq. (3.10) and (3.14) are the fundamental relations governing the eigenvalues of G(y , 9). In the following, 
we explore other properties of these eigenvalues. 

3.2. Characteristic equation of G(y,9) 

By using Eqs. (2.74)-(2.76) and (3.3), one has 


G(v,6) = 

( 2 — cos 9 — in sin 6 2j/(cos 9 — 1) + i(l + v 2 ) sin 9 

2z/(l — cos 9) + i sin 9 (2v 2 — 1) cos 9 — 2v 2 + iv sin 8 

V 3(cos0 — 1) 3v(\ — cos 8) — 3isin0 

— OO < V < +00; — 7T < 8 < 7T 


2(1 + 2v 2 ) 


_ 2iv(2 + u 2 ) a . n ^ 
3 


(1 — cos 9) — — — sin 8 


2u(l + 2v 2 ) 2i(l — v 2 ) . 

— A— -(1 - cos 8) H sin 9 

O O 

2(1 + is 2 ) cos 9 — 1 — 2v 2 + 2 wsin 9 ) 


(3.15) 

(3.16) 

(3.17) 

(3.18) 

The reader may be surprised by the simple result Eq. (3.16). However, by using Eq. (3.3) and the fact 
that each of Qoiy), Q+(y), and Q-iy) has the form dc 1 with c and d being 3x1 column vectors, an 
application of the fundamental definition of determinant (in which the Levi-Civita antisymmetric symbol is 
used) leads to the conclusion that det[G(^, 9)] must be independent of 9 , i.e., clet[G(iy 0)] = det[G(iy 0)]. As 
such, Eq. (3.16) now follows from the fact that G(y, 0) = U (see Eqs. (3.15) and (2.32)) and det(U) = — 1 . 
Hereafter, for simplicity, the arguments v and 9 may be omitted if no confusion would arise. 

Because ay, ay, and 0-3 are the eigenvalues of G, Eq. (3.17) implies that 


It follows from Eq. (3.15) that (i) 

det [G(iy 0)] = — 1 

and (ii) any eigenvalue a of G(u, 9) must be a root of the characteristic equation: 

det [<j I — G(v, 0)] = a 3 + h{ v, 9)<r 2 + h (v, 9)a +1 = 0 

where 

h(v, 9) = f — 1 + Au 2 (l — cos 9) — 2i^sin 9 


<J 3 + ha 2 + ha + 1 = (a — oy )(a - cr 2 )(cr — 03) 
for any complex variable o\ On the other hand, because Eq. (3.10) <t=> 

{01 (iy,9), a 2 (ty,9), a 3 (iy,8)j = ^ 


oy(+ 9) a i (u,9) 0 - 3 ( 1 /, 6 >) 

1/aT, 1/ag, and 1 /ay must also be the eigenvalues. Thus 


a d + her + ha + 1 = [a [a [a 


0-1 


1 


o' 2 


1 


+3 


(3.19) 


(3.20) 


(3.21) 
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for any complex variable a. In the following, Eq. (3.21) will be derived directly from Eq. (3.19) without 
using any other assumption. 

Proof. Eqs. (3.19) 


o 1 cr 2 cr 3 = —1, (7 1 cr 2 + cr 2 tJ 3 + cr 3 fJ i = h, and cr i + cr 2 + cr 3 = —h 


(3.22) 


Eq. (3.9) follows from the relation o\o 2 o^ = —1. Also 010203 = — 1 <t=> Eq. (3.21) is valid if o = 0. 
Let (7^0. Then, by replacing a with l/o in Eq. (3.19), one has 


1 


h 


1 


= 3+ = 2+ = +l— l = -(7i -- <7 2 =-CT 3 

a o o \o 


o <7 

Also, by using the relation oqoqoq = — 1, one has 

a 3 = (— u/<7i)(— ct/<7 2 )(— ct/ct 3 ) 


(3.23) 


(3.24) 


Because the product of the expressions on the left sides of Eq. (3.23) and (3.24) equals to that on the right 
sides, we have 

' 1W 1W ^ (3.25) 


o 3 + ho z + ho + 1 — I o — ) [ o — ) (cr — 




0 1 


1 


02 


1 


(73 


Eq. (3.21) is the complex conjugate form of Eq. (3.25). QED. 
Moreover, according to Eq. (3.18), 


h(—u, 9) = h(v , —9) = h{v, 9), — oo <v< +oo; —7 r < 9 < it 


(3.26) 


Thus Eq. (3.14) can also be derived directly from Eq. (3.19). 

In this section, we will prove the following proposition: 

Proposition 1. \ot(is, 9)\ = 1 for all f and 9, t = 1, 2, 3, and — 7r < 9 < tt, if and only if \u\ < 1/2. 
Proposition 1 can be divided into two parts, i.e., 

Proposition 1(a). \oe(v, 0)| = 1 for all l and 9, i = 1,2,3, and — 7r <9 < 7r, if \v\ < 1/2. 
and 

Proposition 1(b). For any v with \u\ > 1/2, there is a pair of f 0 and 9 0 such that 


(, 0 = 1,2,3, -7T < 9 a < 7r, and \oe o (v,0 o )\ 1 


(3.27) 


A simple proof for Proposition 1(a) will be given in Sec. 3.3. Based on more exhausted developments, another 
proof for Proposition 1(a) and a proof for Proposition 1(b) will be given in Sec. 3.7. 

3.3. A proof for Proposition 1(a) 

First we introduce the following well-established algebraic theorem: 

Theorem 4. Let <7i, cr 2 , . . . , opj> be the distinct roots of the Ath-order algebraic equation 

o N + aicr^ -1 + a 2 o N ~ 2 + . . . + ajv-iu + a,/v = 0 (3.28) 

where aq, a 2 , . . . , ajv are complex constant coefficients and cr is a complex variable. For each £ = 1,2,..., N', 
let me > 1 denote the multiplicity of the root ay. Then 

N' N' 

o N + a\o N1 + a 2 o N 2 + . . . + a^v— icr + ajv = J^[(c r — oe) me and ^ me = N (3.29) 

e=i i=\ 
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According to the above theorem, for any given the roots of the cubic equation Eq. (3.17) must 

fall into one of the following three mutually exclusive cases: (a) there is one triple root (multiplicity = 3); 
(b) there are one double root (multiplicity = 2) and one simple root (multiplicity = 1) and (c) there are 
three simple roots. 

Consider case (a). Then cy = a 2 = 03. Let a Q denote the common value of <Ji, 172, and <73. Then 
Eqs. (3.9) and (3.20) imply that (i) a 0 y£ 0 and (ii) l/oo must also be a triple root of Eq. (3.17). Thus the 
only choice that will not contradict Theorem 4 is that a 0 = 1 /oo, he., |er 0 | = |<7i| = | < 7 - 2 1 = |cr 3 1 = 1. 

Consider case (b). Without any loss of generality, one can assume ey = er 2 y£ cr 3 . Again let a 0 denote 
the common value of ay and 02 ■ Then Eqs. (3.9) and (3.20) imply that (i) cr 0 y£ 0; (ii) cr 3 y£ 0; and (iii) 1/oy 
and 1/oy must also be a double root and a simple root of Eq. (3.17), respectively. Thus the only choice that 
will not contradict Theorem 4 is that <r 0 = 1/oy and cr 3 = l/ay, i.e., |cr 0 1 = |<7i| = |<72 1 = |<r 3 | = 1. 

The conclusions reached above imply the following lemma: 

Lemma 1. For any given ( v,9 ), the roots of Eq. (3.17) must all be of unit magnitude if any one of 


them is a multiple root. 

Thus, to prove Proposition 1(a), we need only to consider case (c), i.e., the case with 

cri y£ cr 2 , oy y£ <r 3 , and cr 2 y£ cr 3 (3.30) 

To proceed, each oy(i/, 9) is expressed in its polar form, i.e., 

ay (v, 9) = re(p, 9)e l ^ e ^' e \ £=1,2,3; — 00 < v < + 00 ; —tt<9<it (3.31) 

where, because of Eq. (3.9) 

dpf 

ry(iy 9) = \ae(y, 9)\ > 0, £=1,2, 3; —00 < v < + 00 ; — 7 r <9 < n (3.32) 

Moreover, for each 9), the corresponding phase angle 4>e(v, 9) is uniquely defined by Eq. (3.31) and 

— 7 r < 4>i(v, 9) < 7r, £ = 1, 2, 3; —00 <v < + 00 ; — 7r < 9 < n (3.33) 

Hereafter, the arguments v and 9 may be dropped from ri(v 7 9) and <fie(i',9) if no confusion would arise. It 
follows from Eqs. (3.31) and (3.32) that 

l/o?=(l/r/)e^, £=1,2,3 (3.34) 

Also, by using Eqs. (3.31)-(3.33), Eqs. (3.30) can be expressed as the following ordered pair inequalities: 

(nAi) ± (n,(j>i) 7 ^ (r3,h), and (r 2 , <£ 2 ) ^ (^3, fo) (3.35) 


The distribution of </>i, </> 2 , and (/> 3 must fall into one of the following mutually exclusive cases: (cl) all 
have distinct values; (c 2 ) two of them have the same value while the third assumes a different value; and 
(c3) all have the same values. In the following, these sub-cases will be discussed separately. 

Consider case (cl) where 

<Ai 4>2, 0 1 7 ^ 03, and 02 y£ 0 3 (3.36) 

Because Eqs. (3.20) and (3.34) imply that (l/r^e*^, £ = 1,2,3, must also be roots of Eq. (3.17), Eq. (3.36) 
implies that the only choice that will not contradict Theorem 4 is that re = 1 /re, i.e., re = 1, £ = 1,2,3. 
Thus, for case (cl), again we have |cti| = |fT 2 1 = |cr 3 | = 1. 

Consider case (c2) where, without any loss of generality, one can assume that 

01 = 02 03 (3.37) 
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Because of (3.35), Eq. (3.37) implies that 

n ± r 2 (3.38) 

By using Eqs. (3.37) and (3.38) along with the fact that (l/r^e*^, £=1,2, 3, must also be roots of Eq. (3.17), 
one concludes that the only choice that will not contradict Theorem 4 is that rir 2 = 1, r\ ^ 1, r 2 1, and 
r 3 = 1. Thus, for case (c2), (i) one of the roots is of unit magnitude while the other two are not; and (ii) 
the product of the magnitudes of the two roots which are not of unit magnitude is one. 

Consider case (c3) where 

<£i = = 03 (3.39) 

Because of Eq. (3.35), Eq. (3.39) implies that 

n ^ r 2 , ri ^ r 3 , and r 2 ^ r 3 (3.40) 


By using an argument similar to that invoked in the discussion of case (c2), one concludes that, for case 
(c3), again (i) one of the roots is of unit magnitude while the other two are not; and (ii) the product of the 
magnitudes of the two roots which are not of unit magnitude is one. 

As a result of the above discussions, we have the following lemma: 

Lemma 2. For any given (i/, 0), the case with at least one of the roots of Eq. (3.17) not being of unit 
magnitude may occur only if it meets the following conditions: (i) one and only one of ri, r 2 , and r 3 is of 
unit magnitude; and (ii) the two roots that are not of unit magnitude share the same phase angle and the 
product of their magnitudes is one. 

Consider any case that meets the conditions referred to in Lemma 2. Then, without any loss of generality, 
one may assume that 

rir 2 = 1, n ^ 1, r 3 = 1, and 0 X = 0 2 = 0 (3.41) 


where 0 denotes the common value of 0i and 0 2 . Moreover, by using Eq. (3.31), Eqs. (3.18) and (3.22) 
imply that 

rir 2 r 3 e iWl+02+03) = -1 (3.42) 

nr 2 e % ^ 1+ ^ + rir 3 e l ^ 1+< ^ >3 ' > + r 2 r 3 e l ^ 2+ ^ = — 1 + 4j/ 2 (1 — cos 9) + 2iusm 9 (3.43) 

and 


rie*^ 1 + r 2 e l ^ 2 + r 3 e l< ^ 3 = 1 — 4^ 2 (1 — cos0) + 2ii/sin0 
Because of Eq. (3.32), Eq. (3.42) 


and 


rir 2 r 3 = 1 


e i(01+02+<i>3) __ 


By using Eqs. (3.45) and (3.46), Eq.(3.43) 


(3.44) 

(3.45) 

(3.46) 


— e*^ 1 H e 1 ^ 2 H e*^ 3 = 1 — 4^ 2 (1 — cos 0) + 2ii/sin 9 (3.47) 

ri r 2 r 3 


Thus Eqs. (3.44)-(3.47) represent all the independent constraints imposed on ry and <f>g, i = 1,2,3. 

Note that: (i) Eq. (3.41) implies Eq. (3.45); and (ii) Eqs. (3.41) and (3.46) imply that 

e i<h = e i<t> 2 = e i<t> and e i0s = ^ e -i(0i+02) = _ e -2i0 (3.48) 


Then, with the aid of Eqs. (3.41) and (3.48), both Eqs. (3.44) and (3.47) reduce to 

f(p)e l<l> — e~ 2 *^ = 1 — 4^ 2 (1 — cos 9) + 2wsin0, — 7r < 0 < 7r; p > 0 and p ^ 1 


(3.49) 

(3.50) 
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where 

f{p) = f P+ P> 0 and p 1 (3.51) 

P 

Eq. (3.50) 

/(p) cos (j) — cos (2^>) = 1 — 4^ 2 (1 — cos0), — n < <fi < tt ; p > 0 and p ^ 1 (3.52) 

and 

/(p) sin 4 + sin(2</>) = 2^sin0, —n < 4 < 7r; p > 0 and p ^ 1 (3.53) 

Thus, given any ( is, 9 ), Eqs. (3.52) and (3.53) must admit a solution for p and </> in the specified domain if 


the case Eq. (3.41) indeed exists. 

To explore Eqs. (3.52) and (3.53), note that (i) 

[/(p) cos 4 > — cos(20)]~ + [/(p) sin</> + sin(20)] 2 = [l — 4 i^ 2 (1 — cos#)] 2 + [2^sin0] 2 (3.54) 

is a direct result of Eqs. (3.52) and (3.53); (ii) 

[f(p) cos 4> - cos(20)] 2 + [/(p) sin 4> + sin(20)] 2 = [/(p) - l] 2 + 2/(p) [1 - cos(3</>)] (3.55) 

and (iii) 

[1 - 4^ 2 (1 - cos 0)] 2 + [2^ sin 6] 2 = 1 - 4^ 2 (1 - 4^ 2 )(1 - cos6») 2 (3.56) 

Next, because (i) the minimum of /(p) in the domain p > 0 occurs at p = 1 and (ii) /( 1) = 2, we have 

/(p) >2 if p > 0 and p yf 1 (3.57) 

Combining Eqs. (3.55) and (3.57), and using the fact that 1 — cos(3^>) > 0 for all />, one has 

[/(p) cos 4> — cos(2</>)]“ + [/(p) sin 4 + sin(2</>)] 2 >1 if p > 0 and p ^ 1 (3.58) 

On the other hand, because 1 — 4^ 2 > 0 if \v\ < 1/2, Eq. (3.56) implies that 

[l — 4^ 2 (1 — cos0)]“ + [2^sin0] 2 < 1 if \v\<l/2 (3.59) 

Combining Eqs. (3.58) and (3.59), one arrives at the conclusion that, for all 0, 

[/(p) cos 4 ~ cos(2^)]“ + [/(p) sin 4+ sin(2^)] 2 > [l — 4 i/ 2 (1 — cos0)] 2 + [2i/sin0]“ (3.60) 


i.e., Eq. (3.54) cannot be satisfied, if (i) \u\ < 1/2; and (ii) p > 0 and p / 1. Because Eq. (3.54) is a direct 
result of Eqs. (3.52) and (3.53), this implies that, for any 9 , Eqs. (3.52) and (3.53) admit no solution for p 
and 4 i n the specified domain, i.e., the case Eq. (3.41) does not exist, if \i/\ < 1/2. In turn, this implies that, 
for all 9 , the roots of Eq. (3.17) are all of unit magnitude if \v\ < 1/2, i.e., Proposition 1(a) has been proved. 
QED. 

Note that, for a reason to be given in Sec. 3.9, by itself Proposition 1(a) does not imply that a( 3) 
scheme is stable when \v\ < 1/2. Next, as a preliminary for later developments, several special cases will be 
discussed in Secs. 3.4 and 3.5. 

3.4. The \v\ = 1/2 case 

Let v = 1/2. Then Eq. (3.17) reduces to 

a 3 — e l6 a 2 — e~ l6 a + 1 = (er — e ld ) (cr — [a + e~ l6 ^ 2 ^j = 0, — tt <9 <tt (y = 1/2) (3.61) 
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Thus the roots of Eq. (3.17) are 

a = <Jo(0) = f e lB and cr = a±{9) A = ±e ~ l8 ^ 2 , —n < 9 < tt {v = 1/2) (3.62) 

On the other hand, by using Eqs. (3.14) and (3.62), one concludes that the roots for the case v = —1/2 are 
a = ao(9) = e ~ l8 and a = cr±{9) = ±e l8 ^ 2 , — n < 9 < tt (^=—1/2) (3.63) 

For each of the above two cases, Eqs. (3.62) and (3.63) imply that the these roots are distinct if 

9^ 0 and \9\ ± — (\v\ = 1/2) (3.64) 

In fact, there are one double root and one simple root if 9 = 0 or \9\ = 2tt/3. Also, because the analytical 
amplification factor is e~ lv8 for any (v,9) (see p.4 of [61]), for the case v — 1/2 ( v — —1/2), one of the roots 
of Eq. (3.17), i.e., o + (60 ( cr +(^)) ) is identical to the analytical amplification factor. 

Consider the plane wave solution 


1 

R" 

'R) 

II 

HO 

{ka ^ 0) 

(3.65) 

The period associated with this solution is 

27T 

T = 

|fca| 


(3.66) 


Let n, the number of total marching steps, and At be chosen such that 


nAt = NT, 


N= 1,2,3,... 


Then, by using Eq. (3.66), one has 


where 


n = 


2ttN 

MR 


9 = kAx ^ 0 

is the variation of phase angle over the interval ax. For the case |^| = 1/2, Eq. (3.68) reduces to 

47T./V 


\ 0 \ 


{nAt = NT ; \v\ = 1/2)) 


Eqs. (3.62), (3.63), and (3.70) imply that 


[a ± {9)] n = cr±{9) =(±1 y 


and 


Thus 


= 1 


= [*±( 6 ) 

On the other hand, Eq. (3.68) implies that 


M0)] n = ko(^)J 

= 1, if n is even 


{e~ ive ) n = 1 


(3.67) 

(3.68) 

(3.69) 

(3.70) 

(3.71) 

(3.72) 

(3.73) 

(3.74) 
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By using an analytical procedure similar to that used in Sec. 5 of [1], one can show that Eqs. (3.72)-(3.74) 
and the fact that e~ w6 is the analytical amplification factor lead to the conclusion that, for the case \u\ = 1/2, 
the numerical solution generated by the a(3) scheme in a simulation involving a periodic boundary condition, 
aside from round-off errors, should be identical to the exact solution if (i) n and At are chosen according to 
Eq. (3.67), and n is even; and (ii) the phase angles of the Fourier components involved in the simulation 
observe the condition Eq. (3.64) (i.e., the three eigenvalues associated with each Fourier component are 
distinct). This prediction has been verified numerically (see Sec. 4). 

Next, a brief discussion on the roots of Eq. (3.17) for the three special cases: (a) v = 0; (b) 9 = 0; and 
(c) 9 = 7T will be given in Sec. 3.5. 

3.5. Three other special cases 

Let v = 0 or 9 = 0. Then Eqs. (3.17) and (3.18) imply that 

CT 3 -cr 2 - cr + 1 = (cr- 1 ) 2 (ct + 1) =0 (3.75) 

i.e., the roots of Eq. (3.17) are 1, 1, and — 1. According to Eqs. (3.31)-(3.33), (i) rq = r 2 = r 3 = 1, and (ii) 
one can assume that <pi = <j) 2 = 0 and <f > 3 = 7 r. 

Let 9 = 7T. Then Eqs. (3.17) and (3.18) imply that 

cr 3 + (8z/ 2 - l)cr 2 + (8z/ 2 - 1)(7 + 1 = (a + 1) [a 2 + 2(4i/ 2 - 1)ct + l] = 0 (3.76) 

i.e., the roots of Eq. (3.17) are a = — 1 (i.e., re = 1 and <j>e = 7r for a value of £), and 

cr = l — 4z/ 2 dz \J&v 2 (2v 2 - 1) (3.77) 


We have 


1 - Av 2 - \j8v 2 (2v 2 - 1) < 1 - 2 = -1 if 2v 2 > 1 


i.e., 


1 - Av 2 - yJSv 2 {2v 2 - 1) 


> 1 if 2v 2 >1 


(3.78) 


Thus the magnitude of at least one root of Eq. (3.17) is greater than one if 9 = 7r and \v\ > l/V^- 
On the other hand, 


1 - 4jA ± \J‘&v 2 (2v 2 - 1) 


1 - Av 2 ± i 2V2 H \/l — 2v 2 


= v / (1-4jA)2 + 8^2(1 _2jA) = 1 if 2v 2 < 1 


(3.79) 


Thus, for the case 9 = tt and \u\ < l/\/2, rq = f2 = r3 = 1. Moreover, for the special case 

9 = 7r and \v\ = \/V2 (3.80) 

Eq. (3.76) reduces to (ct + 1) 3 = 0, i.e., —1 is the triple root of Eq. (3.17). In fact, it will be shown in Sec. 3.6 
that the only possible triple root of unit magnitude for Eq. (3.17) is —1 and it has this root only for the case 
Eq. (3.80). 

Eq. (3.17) has three roots for any given (n,9). In Sec. 3.6, we derive a set of equations governing the 
phase angles of these roots when they all are of unit magnitude. To pave the way, let 


and 


T = {(^, 0)| — 00 < v < +00; — 7T < 9 < tt, and ri(^, 9) = 9) = ^(z^, 9) = 1} 

d'o d = {(v,0) \v ^ 0, 0 < \9\ < 7 r, and r x (i/,0) = r 2 (v,0) = r 3 (z/,0) = 1} 


(3.81) 

(3.82) 
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According to Proposition 1(a) (which has been proved), (y, 9) £ SS? if \v\ < 1/2 and — ir < 9 < n. 

3.6. Phase angle equations for (y, 9) € 'if 

With the aid of Eq. (3.22), and (3.31)-(3.33), one concludes that the condition 

ri{v,9) = r 2 {v,9) = r 3 (v,9) = 1 (3.83) 


<=> the real phase angles <pn{v. 9), l = 1,2,3, satisfy 

g i(01+02+</>3) _ 

e i(0i+0 2 ) _|_ e *(0i+</>3) _|_ e *(02+03) _ _|_ 4jy 2 (i — cos 9) + 2ii/sin0, 

and 


(3.84) 

—oo <v< +oo; —7 r < 9 < tt (3.85) 


e *0i _|_ e *02 _|_ e *03 _ i _ 4^2(1 _ cos 9) + 2ivsin6, — oo <v< +oo; — 7r < 9 < n (3.86) 

Because (i) Eq. (3.84) implies that 

gi(01+0 2 ) _ _g-#3 g*(01+03) __ _g-*02 g*(02+03) _ _g-*01 ^3 

and (ii) the complex conjugate of Eq. (3.86) is 

g-*0i _|_ e -*^2 _|_ e -*03 _ 4 _ 4 j, 2 (i _ cos 9) — 2iusin9 (3.88) 

one can see easily that Eq. (3.85) is a result of Eqs. (3.84) and (3.86). Thus Eq. (3.83) Eqs. (3.84) and 
(3.86). 

At this juncture, we will prove that the only possible triple root of unit magnitude for Eq. (3.17) is — 1 
and it has this root only for the case Eq. (3.80). 

Proof. Let Eq. (3.17) have a triple root of unit magnitude. Then (;/, 9) £ T and fa = fa — fa- With 
the aid of Eqs. (3.33), in turn Eq. (3.84) implies that either (a) 

fa = fa = fa = tt (3.89) 

or (b) 

fa = fa = fa = ±tt/3 (3.90) 

For case (a), by substituting Eq. (3.89) into Eq. (3.86), one has 

is 2 (l — cos 9) = 1 and usin0 = 0 (3.91) 

which <=> cos9 = —1 and \v\ = 1 / a/ 2- Because — 7r < 9 < w, in turn Eq. (3.91) Eq. (3.80), i.e., case (a) is 
the same case defined by Eq. (3.80). 

For case (b), by Substituting Eq. (3.90) into Eq. (3.86), one has 

(3/2) (1 ± \/3 i) = 1 — 4u 2 (1 — cos 9) + 2ivsm0 (3.92) 


By taking the real part of Eq. (3.92), one arrives at the result 

v 2 (l — cos0) = —1/8 (3.93) 

Because v 2 (l — cos 9) > 0 for any fa, 9) and —1/8 < 0, Eq. (3.93) cannot be true and therefore case (b) does 
not exist. QED. 

Let (y,9) € T. Then Eqs. (3.84) and (3.86) are valid. By eliminating e *^ 3 from Eqs. (3.84) and (3.86), 
one has 

e *<A _|_ e *02 _ g-*(0i+02) _ i _ 4^(1 _ cos 9) + 2ivsin9 (3.94) 
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By separating the real and imaginary parts of Eq. (3.94), we have 

cos 0 i + cos 0 2 — cos(0i + 0 2 ) — 1 = — 4^ 2 (1 — cos 9) (3.95) 

and 

sin 0 1 + sin 02 + sin(0i + 0 2 ) = 2r> sin 9 (3.96) 

Similarly, one can show that 

cos 0 2 + cos 4>3 — cos (</>2 + <^ 3 ) — 1 = — 4^ 2 (1 — cos 9) (3.97) 

sin 0 2 + sin 0 3 + sin (02 + $ 3 ) = 2v sin 9 (3.98) 

cos 0 3 + cos 0i — cos ((/>3 + 0i) — 1 = — 4^ 2 (1 — cos 9) (3.99) 

and 

sin 03 + sin + sin(0 3 + 0i) = 2^ sin 9 (3.100) 

At this juncture, note that Eqs. (3.95), (3.97), and (3.99) imply that 

^ 2 (1 — cos 9) = 0 if fie = 0 for any £ = 1,2,3 (3.101) 

Because — 7 r < 9 < n, Eq. (3.101) implies that at least one of the two cases (a) v = 0 and (b) 9 = 0 must 
occur if cj>e = 0 for any £=1,2,3. It is shown in Sec. 3.5 that, for v = 0 or 9 = 0, indeed <f>e = 0 for two 
different values of l. 

On the other hand, Eqs. (3.96), (3.98), and (3.100) imply that 

^sin0 = 0 if <f>e = 7 r for any l = 1, 2, 3 (3.102) 

Because — 7 r < 9 < it, Eq. (3.102) implies that at least one of the three cases (a) v = 0, (b) 9 = 0, (c) 9 = n 
must occur if 0^ = 7 r for any £=1,2, 3. It is also shown in Sec. 3.5 that, for any of above three cases, indeed 
= 7 r for a value of £. 

Using the above results along with Eqs. (3.33) and (3.82), one arrives at the important conclusion that 

0 < |0*(i/,0)| < 7T, £=1,2,3, if (i/,0)etf o (3.103) 

To eliminate 02 , let (i) Eq. (3.95) be multiplied by (l+cos0i); and (ii) Eq. (3.96) be multiplied by sin0i. 
After subtracting the resulting equations from each other, a rearrangement using elementary trigonometry 
yields 

sin 2 0i — 2 j/ 2 (1 — cos0)(l + cos0i) — ^sin0sin0i = 0 (3.104) 

By applying similar manipulations over Eqs. (3.97), (3.98), (3.99), and (3.100), one can show that Eq. (3.104) 
remains valid if 0i is replaced by 0 2 or 0 3 . Thus, for any (^,0) € db we have 

F(v,0,<j> e )= 0, £=1,2,3 

where 

F{y , 9, 0) = f sin 2 0 — 2 iA( 1 — cos 0)(1 + cos 0) — v sin 9 sin 0 
— OO < V < +OO; — 7T < 9 < 7T; — 7T < 0 < 7T 
Eq. (3.106) implies that, for all (v,9) with —00 < v < +00 and —ir < 9 < ir, we have 

F{-v, 9, 0) = F{v, - 9 , 0) = F(v, 9 , -0) 

F(±l/2, 9, ±9) = F(± 1/2, 9 , =f 6»/2) = F(±l/2, 0, ir =F 9/2) = 0 


(3.105) 

(3.106) 

(3.107) 

(3.108) 
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and 


F(v, 9 ,tt) =0 (3.109) 

Eqs. (3.105) and (3.107) are consistent with Eqs. (3.14) and (3.31) while Eq. (3.108) is consistent with the 
special results Eqs. (3.62) and (3.63). On the other hand, because (i) sin7r = l+cos7r = 0, and (ii) Eq. (3.104) 
is obtained from a subtraction of two expressions which result from multiplying Eqs. (3.95) and (3.96) with 
(1+cos 0i) and sin 0 i, respectively, the fact that Eq. (3.109) is true for all (i/, 9 ), — oo < v < +oo; — tt < 9 < tt. 
is an artificial result accidently introduced in the derivation of Eq. (3.105). 

According to the above discussions, given any (v,9) € 4', the phase angle 0 of any root of Eq. (3.17) 
must satisfy 

F(v , 9 , 0) = 0, —oo <v< +oo; — tt < 9 <tt\ —tt < 0 < tt (3.110) 

Recall that the analytical amplification factor is given by e~ w6 . Thus it is expected that 0 = — ^0 should 
be a good approximated solution to Eq. (3.110) when \6\ is small (i.e., when the solution Eq. (3.65) has a 
very small variation over the spatial interval Ax and thus it is closely approximated by a discrete solution). 
In fact, with the aid of Eq. (3.106) and the Taylor’s expansions 


x 


x 


x ’ 


x 


s inz = z-— + — + 0(x 9 ) and cosz = 1 + 0(a: 8 ) 


6 120 5040 


24 720 


(3.111) 


one has 


„ „ w2$6 

F{v 1 9 1 —v9) = sin 2 (v9) — 2^ 2 (1 — cos0)[l + cos(^0)] + z'sin0sin(z'0) = (Av 2 — 1)(^ 2 — 1) + O(0 8 ) (3.112) 

oou 

Because Av 2 — 1 = 0 <t=> \v\ = 1/2, the above result is consistent with the fact that F(±l/2, 9, =p0/2) = 0, 
which was presented as part of Eq. (3.108). 

Given any ( v,9 ), Newton’s iterative procedure for obtaining a root 4> of Eq. (3.110) is defined by 


^n+l = 


F(v,9,^) 

F<f>{v, 9, <t> n ) ’ 


n = 0, 1, 2, 3, . . . 


(3.113) 


where (i) <j) n is the ?ith iterative value of (f> and (ii) 

FAv, 9, (j)) = f = sin(2t/>) + 2v 2 (1 — cos 9) sin (j> — v sin 9 cos 4> (3.114) 

ocp 

For a given (y,9) £ T, the phase angle <j> of any crt(y,9) must satisfy Eq. (3.110). Moreover, according 
to Eq. (3.103), 0 < \<f>\ < tt if v yf 0 and 0 < \9\ < tt. As a preliminary to the proof to be given in Sec. 3.7, 
the equation F(y, 9, <j>) = 0 will be cast into a cubic equation assuming 

v ^ 0 and 0 < \9\, \<j>\ < ir (3.115) 

As a result of Eq. (3.115), we have 


1 + cos (j)^ 0 and v sin 9 yf 0 
With the aid of Eqs. (3.106) and (3.116), Eq. (3.110) implies that 
1 — cos (j> — 2v 2 (1 — cos 9) sin (j) 


v sin 9 


1 + cos ( 


= 0, 


^ ± 0; 0 < |0|, \(j)\ < tt 


(3.116) 


(3.117) 


Because 


1 — COS (j) = 


2 tan 2 (0/2) 

1 + tan 2 (0/2) ’ 


sin 0 
1 + cos 0 


= tan(0/2), 


and 


1 — cos 9 
sin 9 


tan(0/2), 0 < \9\, |0| < tt (3.118) 
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Eq. (3.117) <^> 

t 3 + 2 zHan(0/2) t 2 + t + 2z/tan(0/2) = 0, 0 0; 0 < 101 < n (3.119) 

v sm 9 \ 

where r is related to <j> through the one-to-one relation 

t = f tan(c/>/2), 0 < 1 0| < 7r (3.120) 

It has been shown that, for any (v,9) £ 'h 0 , each real phase angle 9), £ = 1,2,3 must satisfy 
Eq. (3.119) through the one-to-one relation r = tan 9)/2}. As such, assuming (i) (u, 9) £ 'ho and (ii) 
£ = 1,2,3 are distinct, Theorem 4 implies that the roots of Eq. (3.119) are also real and distinct, 
and can be indexed such that 

ti(v, 6) = tan [<f>t(u, 9)/2] , £=1,2,3 (3.121) 

However, for the case in which <j>e(v, 9), l = 1,2, 3, are not distinct, two or three of the phase angles that share 
a common value may be linked with a real root of Eq. (3.119) through a relation in the form of Eq. (3.120). 
As such there is a possibility that one or two roots of Eq. (3.119) may not be linked with any 9) in the 
form of Eq. (3.120) or any way whatsoever. In the following, this possibility will be ruled out. In fact, we 
will prove the following proposition: 

Proposition 2. Let (i) v 0 0 and 0 < \6\ < 7r, and (ii) rt{v, 9) and (j>e(v, 9), £=1,2, 3, be defined using 
Eqs. (3.31)-(3.33). Then {y,9) £ 'k 0 , i.e., Eq. (3.83) is true, if and only if the roots of Eq. (3.119) are all 
real. Moreover, these real roots can be indexed such that they and the real phase angles <00, 9), £ = 1, 2, 3, 
are related through Eq. (3.121). 

Proof. Let n{v,9), £ = 1,2,3, denote the roots of Eq. (3.119) where (u, 9) is only subjected to the 
condition v 0 0 and 0 < \8\ < i r. Then, no matter how these roots are assigned the indices £ = 1,2,3, we 


have 

t 3 + 2 j/tan(#/2) r 2 + r + 2zHan(0/2) = (r — n)(r — T2)(t — 73), ^0 0; 0 < \9\ < tt (3.122) 

v sm 9 

Hereafter, for simplicity, the arguments v and 9 may be dropped from th(v, 6) and tan \<!>i{y, 9) /2], £=1,2, 3. 
Eq. (3.122) 

T1T2T3 = — 2zHan(#/2), ^0 0; 0 < \9\ < tt (3.123) 

T1T2 + t 2 t 3 + t 3 t 1 = 1 , v 0 0; 0 < \6\ < tt (3.124) 

and 

ti + t 2 + t 3 = 2 — - — vtan(6/2) , v 0 0; 0 < 101 < ir (3.125) 

v sm 9 

Because vtan(9/2) 0 0 if v ^ 0 and 0 < \9\ < 7 r, Eq. (3.123) implies that 

ti 0 0, r 2 0 0, and T3 0 0, (3.126) 

Moreover, it follows from Eqs. (3.123)-(3.125) that the roots of Eq. (3.119) are all real and can be indexed 
such that they are related to (j>e{v, 9), £=1,2, 3, through Eq. (3.121), if and only if 

tan(<^i/2) tan((/> 2 /2) tan(<fo/2) = —2utan(9/2), ^0 0; 0 < \9\ < 7r (3.127) 

tan((()i/2) tan(02/2) + tan(02/2) tan(03/2) + tan(</>3/2) tan(</>i/2) = 1, ^0 0; 0 < \9\ < 7r (3.128) 
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and 


tan(</>i/2) + tai\(fa/2) + tan(</> 3 /2) = 2 


1 

^ sin 9 


vtan(9/2) , 


yf 0; 0 < |0| < 7T 


(3.129) 


In the following, first we will show that Eq. (3.119) has only real roots and they can be specified by Eq. (3.121), 
i.e., Eqs. (3.127)-(3.129) are true, if ( v,9 ) e 'Iq,. 

As a preliminary, first note that Eq. (3.83) Eqs. (3.84) and (3.86) (a conclusion reached following 
Eq. (3.88)). Next, because of Eq. (3.33), Eq. (3.84) either (i) Eq. (3.89), or (ii) 


fa + ^2 + fa = ±7T 


(3.130) 


It was shown earlier that Eq. (3.89) can occur only for the case Eq. (3.80). Because this case is ruled out 
by the condition 0 < |0| < 7r, Eq. (3.84) Eq. (3.130) if 0 < |0| < 7r. Moreover, assuming Eq. (3.130), one 
can see easily that Eq. (3.86) <t=> Eq. (3.94) <t=> Eqs. (3.95) and (3.96). Thus one concludes that Eq. (3.83) <t=> 
Eqs. (3.95), (3.96), and (3.130) if 0 < \6\ < 7 r. Note that Eqs. (3.97)-(3.100) are trivial results of Eqs. (3.95), 
(3.96), and (3.130). 

Let fa, 9) G $ 0 . Then according to Eqs. (3.82) and (3.103), and the above discussions, we have (i) 
Eqs. (3.95), and (3.96), and (ii) 


v ^ 0, 0 < \9\ < tv, 0 < \fa\, \fa\, |< 5^3 1 < 7T, and fa + fa + fa = ±7r 

To prove Eq. (3.128), note that the last two conditions given in Eq. (3.131) imply that 

tsm(fa/2) = tan 


(3.131) 


(3.132) 


(Note: t&n(fa/2) is undefined when fa = ±7t,±37t, ±57t, . . .. However these undefined cases are ruled out 
by the condition 0 < | <^>3 1 < 7r.) Eq. (3.128) follows immediately from Eq. (3.132) and the relation 


, /,7T 4>i+4>2\_ ( 4>1 + 4>2\ _ 1 - tan(^i/2)tan((()2/2) 

V 2 2 )- V 2 /” tan(^/2)+tanfe/2) 


(3.133) 


(Note: Because of the last two conditions given in Eq. (3.131), all terms and expressions which appear in 
Eq. (3.133) is well defined.) 

To prove Eqs. (3.127) and (3.129), note that, because of Eq. (3.131), the term 2^sin 9 on the right side 
of Eq. (3.96) is nonzero. Thus the expression on the left side is also nonzero, i.e., 


sin fa + sin fa + s\n(fa + fa) yf 0 


(3.134) 


As such Eq. (3.96) <t=> 


sin fa + sin fa + sin(0i + fa) v sin 6 


(3.135) 


if Eq. (3.131) is assumed. Also by dividing Eq. (3.95) over (3.96), and using the last identity presented in 
Eq. (3.118), one has 

cos fa + cos fa — cos(fa + fa) — 1 


sin fa + sin fa + s\n(fa + fa) 
Adding Eq. (3.135) to Eq. (3.136), we have 

cos fa + cos fa — cos(fa + fa) + 3 


sin fa + sin fa + sin(0i + fa) 


= 2 


= — 2^tan(0/2) 


— iztan(0/2) 


v sin ( 


(3.136) 


(3.137) 
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Furthermore, by assuming Eqs. (3.130) and (3.134), it is shown in Appendix B that 


cos 0 i + cos 0 2 — cos(0i + 0 2 ) — 1 
sin (f>i + sin 0 2 + sin(0! + 0 2 ) 


tan(0i/2) tan(0 2 /2) tan(03/2) 


(3.138) 


and 


COS 01 + COS 02 — cos(0i + 02 ) + 3 
sin 0! + sin 0 2 + sin(0! + 0 2 ) 


tan(0i/2) + tan(0 2 /2) + tan(03/2) 


(3.139) 


Eqs. (3.127) and (3.129) now follow from Eqs. (3.136)-(3.139). Thus we have shown that the roots of 
Eq. (3.119) are the real roots specified by Eq. (3.121), if (i/, 0)'F O . 

Next consider any (i/, 0) such that (i) v ^ 0 and 0 < |0| < 7 r; and (ii) the roots of Eq. (3.119) are all 
real. In the following we will complete the proof by showing that, for such a ( v , 0), (i) Eq. (3.83) is true, i.e., 
(v, 9) £ T 0 ; and (ii) the real roots of Eq. (3.119) can be indexed such that they and the real phase angles 
0^(z', 9), £ = 1,2, 3, are related by Eq. (3.121). 

Let the real roots be denoted by Te(u,9), £ = 1,2,3. Then we have Eqs. (3.123)-(3.126). Because of 
Eq. (3.126), for each iy(iy 0), there exists one and only one p^, 0) such that 


0 < |pi(i2,0)|, |P2(z',0)|, |P 3 iy-> 0 ) | < 7 T 


and 


n(v,9) = tan[p^(z/,0)/2] , ^ = 1,2,3 

Substituting Eq. (3.141) into Eq. (3.123)-(3.125) and dropping the arguments v and 0, we have 
tan(pi/2) tan(p 2 /2) tan(p 3 / 2 ) = — 2j/tan(0/2), v ^ 0; 0 < |0| < 7r 


(3.140) 

(3.141) 

(3.142) 


tan(pi/2) tan(p 2 /2) + tan(p 2 /2) tan(p3/2) + tan(p3/2) tan(pi/2) = 1, v ^ 0; 0 < |0| < 7r (3.143) 


and 


tan(p!/2) + tan(p 2 /2) + tan(p 3 / 2 ) = 2 


1 


v sin 9 


— ^tan(0/2) 


v ^ 0; 0 < |0| < 7T 


(3.144) 


Because of Eq. (3.140), at least two of p 1 , p 2 , and P 3 must both be positive or negative. Let pi and 
P 2 be both positive or negative, i.e., P 1 P 2 > 0. Then Eq. (3.140) also implies that tan(pi/2) tan(p 2 /2) > 0. 
In turn, we have 

pi+p 2 ^0 (3.145) 


and 


tan(pi/2) + tan(p 2 /2) ^ 0 


(3.146) 


(Note: Recall that (a + b) 2 = a 2 + b 2 + 2ab. Thus (a + b) 2 > 0, i.e., a + b ^ 0, if ab > 0.). Next, by combining 
Eqs. (3.140) and (3.145), one has 

I Pl + P2 


0 < 


< 7T 


and therefore 


. ,Pl + P2, , „ 

sm( ) 0 


By using Eq. (3.140) and (3.148), one can show easily that 

■ Vi + P 2 \ _ 1 - tan(pi/2) tan(p 2 /2) 
tan(pi/2) + tan(p 2 / 2 ) 


(3.147) 

(3.148) 

(3.149) 
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Moreover, with the aid of Eq. (3.146), Eq. (3.143) implies that 


tan(<p 3 /2) = 


1 — tan(<pi/2) tan(<^2/2) 


tan((/?i/2) + tan(t/? 2 /2) 

Combining Eqs. (3.149) and (3.150), one arrives at the conclusion that 

tan(<p 3 /2) = cot 

Eq. (3.151) implies that pi, ip 2 , and y? 3 must satisfy one of the following conditions: 

Pi + P2 + PZ = W7T, TO = ±1, ±3, ±5, . . . 

Because 0 < \pi + p 2 + <p 3 | < 37r is required by Eq. (3.140), Eq. (3.152) now implies that 

Pi + P2 + P3 = ±7T 


(3.150) 


(3.151) 


(3.152) 


(3.153) 


Note that, using similar arguments, we will arrive at the same conclusion Eq. (3.153) if, instead of assuming 
pip 2 > 0 at the beginning, we assume p 2 pz > 0 or pzpi > 0- 

To proceed, note that: (i) by combining Eqs. (3.140) and (3.153) with the assumptions v ^ 0 and 
0 < \6\ < 7 r, we have 

0, 0 < \9\ < IT, 0 < \pi\, \p 2 \, |v?3 1 < 7T ; and pi + p 2 + Pz = ±?r (3.154) 


and (ii) 

is a result of Eq. 
Eq. (3.140), i.e., 


sin</?i + sin ip 2 + sin(y?i + (p 2 ) = ±4cos(y>i/2) cos(p 2 /2) cos(y> 3 /2) (3.155) 

(3.153) (see the proof given in Appendix B). By combining Eq. (3.155) with a result of 

cos(^/2)>0, £=1,2,3 (3.156) 


one concludes that 


sin^i + sin + sin(<pi + p 2 ) ^ 0 (3.157) 

is a result of Eq. (3.140) and (3.153). Using arguments similar to those used in the proof of Eqs. (3.138) and 
(3.139) (see Appendix B), one can prove that 


cos pi + cos p 2 — cos(y>i + p 2 ) — 1 
sin (px + sin p 2 + sin(<pi + p 2 ) 


tan(t/?i/2) tan(^2/2) tan(<y2 3 /2) 


and 


cos pi + cos p 2 — cos(y>i + p 2 ) + 3 
sin <pi + sin p 2 + sin(<^i + p 2 ) 


tan(<pi/2) + tan(^ 2 /2) + tan(y> 3 /2) 


follows from Eqs. (3.153) and (3.157). 

Eqs. (3.142), (3.144), (3.158), and (3.159) now imply that 


and 


c» y ,+co sya -co s to + y2 )-l = _ 2l/tanW2) 
sin (p 1 + sin p 2 + sin(<£i + p 2 ) 


cos pi + cos p 2 — cos(^i + p 2 ) + 3 
sin ipi + sin <p 2 + sin(<^i + <p 2 ) 


2 


1 

v sin 9 


z/tan(#/2) 


(3.158) 


(3.159) 


(3.160) 


(3.161) 
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By subtracting Eq. (3.160) from Eq. (3.161), and then taking the reciprocals of the expressions on the both 
sides of the resulting equation, one has 

sin + sin <^2 + sin(<^i + (^ 2 ) = 2^sin0 (3.162) 

Moreover, by substituting Eq. (3.162) into Eq. (3.160) and using the last identity presented in Eq. (3.118), 
one has 

cos <pi ± cos <p2 — cos(y>i + p 2 ) — 1 = — 4i/ 2 (1 — cos 9) (3.163) 

At this juncture, note that Eqs. (3.153), (3.154), (3.162), (3.163) will become Eqs. (3.130), (3.131), (3.96), 
and (3.95), respectively, if the symbols p\, P2, and p^ in the former equations are replaced by <fii, (fi 2 , and 
c t>3 , respectively. 

Next, by using Eq. (3.153), (3.162), and (3.163), one has 

e i(<Pi+<P2+<p 3 ) = _i (3.164) 

and 

e z(fil ± e l<P2 ± e z(fi3 = 1 — 4 iA (1 — cos0) + 2iusm9 (3.165) 

Because Eq. (3.85) is a result of Eqs. (3.84) and (3.86), one can see that 

e i( Vl+V2 ) + e i( V i+ V3 ) + e i(v 2 +vs) = _i + 4 z ^ 2 ( 1 - cos 9) ± 2ivs\w 9 (3.166) 

is a result of Eqs. (3.164) and (3.165). By comparing Eqs. (3.164)-(3.166) with Eq. (3.22), one concludes 
that the roots of Eq. (3.17) are e lvt , £ = 1, 2, 3. Thus, according to Eqs. (3.31)-(3.33) and (3.140), Eq. (3.83) 

def 

is true, and one can choose (fit = pt, £ = 1,2, 3. As such, it has been shown that, for any ( v , 9) such that (i) 
v 0 and 0 < |0| < n; and (ii) the roots of Eq. (3.119) are all real, Eq. (3.83) is true and, through a proper 
indexing, the real roots Tg(u, 9), £=1,2, 3, of Eq. (3.119) and the real phase angle , 9), £ = 1,2, 3, of the 
roots to Eq. (3.17) are related through Eq. (3.121). Thus the proof for Proposition 2 is completed. QED. 

According to Eqs. (3.62) and (3.63), for the special cases v = ±1/2, (i) Eq. (3.83) is true in the domain 
— 7 r < 9 < 7r; and (ii) in the domain 0 < \9\ < n, the phase angles (fit , € = 1,2,3 (which are subjected to the 
condition Eq. (3.33)), can be chosen as 


(fi! = ±9, cfi2 = T9/2, and </>3 = { ^ + if 0 > 0 > -tt ’ 17 = ±1/2; ° < ^ < * ( 3 ‘ 167 ) 

(Note: Hereafter, for Eq. (3.167) and similar equations associated with the special cases v = ±1/2, each 
equation is valid when the upper (lower) signs are taken uniformly.) According to Eq. (3.121) and (3.167), 
the roots to Eq. (3.119) for the current special cases are 

n = ±tan(0/2), T2 = =Ftan(0/4), and T 3 = ±cot(0/4) v = ±1/2; 0 < |0| < 7 r (3.168) 

In fact, by using the relations 

— tan(6 l /2) ± tan(0/4) — cot(0/4) = tan(0/2) T — -, 0 < \9\ < 7 r (3.169) 

sin 9 

and 

tan(0/2) [cot(0/4) — tan(0/4)] = 2, 0 < |0| < 7 r (3.170) 

(which are proved in Appendix B), Eq. (3.168) implies that 




4 

sin 9 


tan(0/2) , 


T\T 2 ± T 2 T 3 ± T3T1 


1, and T1T2T3 = =Ftan(0/2) 


(3.171) 
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In other words, T£, l = 1,2,3, given in Eq. (3.168) indeed satisfy Eqs. (3.123)-(3.125) for the special cases 

v = ±1/2. 

According to Eq. (3.168), we have (i) 

r 2 ^ T\ and r 2 ^ 73 , v = ±1/2; 0 < |0| < tt (3.172) 

(ii) 

ri ^ 73 if v = ±1/2, 0 < |0| < 7 r, and |0| 27 t/3 (3.173) 

and (iii) 

7"i — T 3 if 1 / = ±1/2 and |0| = 2tt/3 (3.174) 

As a result, for the special cases v = ±1/2, Eq. (3.119) has: (i) three distinct real roots if 0 < |0| < 7r and 
\9\ 27r/3; and (ii) one doubt real root (i.e., iq = 73 ) and one simple real root (i.e., r 2 ) if |0| = 27r/3. 

3.7. Proof for Proposition 1 

As a preliminary, we introduce the following well-known theorem: 

Theorem 5. Consider the cubic equation 

r 3 ± a 2 r 2 ± air ± oq = 0 (3.175) 


where a 0 , ai, and a 2 are real coefficients. Let 


def 



(«2) 2 def «1 02 — 3ao 

9 ’ r " 6 


M 3 

27 ’ 


and D a ^q 3 


(3.176) 


Then Eq. (3.175) has: (i) one real root and a pair of complex conjugate roots if D > 0; (ii) three real roots 
and at least two are equal if D = 0; and (iii) three distinct real roots if D < 0. 

For each ( 1 /, 0), Eq. (3.119) is a special case of Eq. (3.175) with 


a 0 = f 2z/tan(0/2), a\ = f 1, and a 2 = f 2 


z,tan(0/2) — 


1 


v sin 9 

Thus, for each {v,6), the discriminant D associated with Eq. (3.119) has the form 


, j/^ 0; 0< |0| <ir (3.177) 


D(v,9) = \ ^ ^ \v tan(0/2) l — 

3 9 v sin 


i[2i/tan(0/2) ± 

3 L v sin 0 


27 L 


i/tan(0/2) — 


' sin 9 


1 4 

27 _ 27 L 
1 


z^tan(0/2) — 


1 


v sin ( 


2 16 

+ 81L 


z^tan(0/2) — 


1 


v sin 9 


9 L 


2i/tan(0/2) 

1 


1 


27 v 2 sin 2 9 

8v 2 tan 2 (0/2) 
27 


1 + 


1 - 


2 16 

v sin 0 J + 8lL 
16tan(0/2) 


2^tan(0/2) 


v sm ( 


z/tan(0/2) — 


1 l 3 


v sm U J 


sin 0 

6tan(0/2) 


1 

27 


1 + 


48 tan 2 (0/2) 2Otan(0/2) 


sin 2 0 


sm( 


smf 


± ^ - tan 4 (0/2), i/ y^ 0; 0 < |0| < tt 


Let 


def o 
S = 


(3.178) 

(3.179) 
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Then Eq. (3.178) and the relation 


tan (0/2) = sin(0/2)/cos(0/2) = 1 2 , , . 

sin 9 2sin(0/2) cos(0/2) 2 1 7 h 


0 < |0| < 7T 


(3.180) 


imply that 


r](s,0) = f — — — ^ ^ ^ = [tan 4 (0/2) sin 2 0] s 3 ± - tan 2 (0/2) sin 2 0[l — 3 sec 2 (0/2)] s 2 

+ sin 2 0[l + 12sec 4 (0/2) + 10 sec 2 (0/2)1 s — [l + 8 sec 2 (0/2)1 , s > 0; 0 < 101 < n 
16 16 

For the special cases v — ±1/2, i.e., s = 1/4, r](s, 9 ) reduces to 


(3.181) 


77(1/4, 9) = — | tan 4 (9/2) sin 2 9 + 2 tan 2 (0/2) sin 2 0[l — 3 sec 2 (0/2)] 

+ sin 2 0[l ± 12sec 4 (0/2) ± lOsec 2 (0/2)] — 4[l + 8sec 2 (0/2)] j, 0 < |0| < n 
By using trigonometric relations such as 

tan 2 (0/2) = sec 2 (0/2) — 1 and sec 2 (0/2) sin 2 0 = 4 sin 2 (0/2) 


Eq. (3.182) can be simplified as 


0(1/4, 0) 


1 

16 


4 cos 2 (0/2) — 1 
cos(0/2) 


0 < |0| < 7T 


(3.182) 


(3.183) 


(3.184) 


Because (i) |0| < tt <t=> |0/2| < 7r/2; (ii) 4cos 2 (0/2) = 1 <t=> cos(0/2) = 1/2 if |0/2| < 7 t/ 2; and (iii) 
cos(0/2) = 1/2 <t=> |0/2| = 7 t/ 3 if |0/2| < n/2, Eq.(3.184) implies that 


0(1/4 ,0) 


< 0 if 0 < |0| < 7r and \9\ ^ 2tt/3 
= 0 if |0| = 2 tt/3 


(3.185) 


With the aid of Eq. (3.179) and the definition of r/(s,9) given in Eq. (3.181), Eq. (3.185) and Theorem 5 
imply that, for the case with v = ±1/2, Eq. (3.119) has: (i) three distinct real roots if 0 < |0| < 7r and 
|0| 7 ^ 27t/3; and (ii) three real roots and at least two are equal if |0| = 27t/3. This result is consistent with 
the conclusion reached following Eq. (3.174). 

Let 

4(0) = f 3tan 4 (0/2) sin 2 0, 0 < |0| < 7r (3.186) 



B{9) = f tan 2 (0/2) sin 2 0 [l — 3sec 2 (0/2)] , 

0 < 0 < 7T 

(3.187) 

and 

C{9) d = i_sin 2 0[l±12sec 4 (0/2)±lOsec 2 (0/2)] , 

0 < 1 0 1 < 7T 

(3.188) 

Then (i) 

4(0) > 0, B{9) < 0, and C{9) > 0, 

0 < 0 < 7T 

(3.189) 

(ii) 





dy(s, 9) 

ds 


= A(9)s 2 +B(9)s+C(9) = 4(0) 


s± 


B(9) I 2 44(0)C(0) — [B(9)] 2 


24(0) 


4[4(0)] 2 


, s > 0; 0 < |0| < tt (3.190) 
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and (iii) because sec 2 (9/2) > 1 , 0 < |0| < tt, 

4 A(0)C(0) - [B(9)} 2 = tan 4 (0/2) sin 4 0 [(27/2) sec 2 (0/2) - (1/4)] >0, 0 < |0| < tt (3.191) 

Eqs. (3.189)-(3.191) now imply that 

df ](s,e) > (I s > 0; 0 < |0| < 7r (3.192) 

as 

Combining Eqs. (3.185) and (3.192), one arrives at the conclusion that (i) 

r](s, 0) < 0 if 0 < s < 1/4 and 0 < \9\ < 7r (3.193) 

and (ii) 

rj(s,9) >0 if s > 1/4 and \9\ = 27t/3 (3.194) 

By using Eqs. (3.185), (3.193), (3.179), and (3.181), one concludes that 

( < 0 if 0 < \v\ < 1/2 and 0 < \9\ < n 

D{v, 9) l < 0 if \v\ = 1/2, 0 < \9\ < tt, and \9\ ± 2tt/3 (3.195) 

[ = 0 if \v\ = 1/2 and \9\ = 27r/3 


With the aid of Theorem 5, Eq. (3.195) infers that the roots of Eq. (3.119) are real and distinct if (i) 
0 < \v\ < 1/2 and 0 < \9\ < tt, and also if (ii) \u\ = 1/2, 0 < \9\ < n, and \9\ ^ 2ir/3. Moreover, it infers that 
these roots are real and at least two of them are equal if \v\ = 1/2 and \9\ = 27t/ 3. By using Proposition 2, 
in turn, one concludes that <Jt{v,9), l = 1,2,3, are distinct and of unit magnitude if (i) 0 < \v\ < 1/2 and 
0 < |0| < 7 r, and also if (ii) \v\ = 1/2, 0 < \9\ < i r, and \9\ ^ 2tt/3. Moreover, one concludes that cre(i>, 9), 
£= 1 , 2 , 3 , are of unit magnitude and at least two of them are equal if \v\ = 1/2 and \9\ = 2tt/3. 

Proposition 1(a) now follows from the above conclusions and several results obtained in Sec. 3.5, i.e., 

(i) the roots of Eq. (3.17) are 1, 1, and — 1 if v = 0 or 9 = 0; and (ii) the roots of Eq. (3.17) are —1 and 
two distinct complex conjugate numbers of unit magnitude if 0 < \u\ < \f\f2 and 9 = 7r (see Eq. (3.77) and 
(3.79)). On the other hand, Proposition 1(b) is a trivial result of (i) Eqs. (3.194), (3.181), and ( 3 . 179 ); (ii) 
Theorem 5; and (iii) Proposition 2. Thus the proof for Proposition 1 is completed. QED. 

In Sec. 3.8, we introduce and prove Proposition 3, which along with the results presented in Sec. 3.5, 
defines the sets of iy,9) for which Eq. (3.17), respectively, has (i) three distinct roots of unit magnitude, 

(ii) one double root of unit magnitude and one simple root of unit magnitude, (iii) one triple root of unit 
magnitude, and (iv) at least one root not of unit magnitude. 

3.8. Proposition 3 

Note that, for any given 9 with 0 < \9\ < tt, (i) Eqs. (3.192) and (3.193) imply that r](s,9) is a strictly 
monotonically increasing function of s in the domain s > 0 and it becomes negative uniformly in the domain 
0 < s < 1/4; and (ii) the coefficient tan 4 (0/2) sin 2 9 of the third-order term in s in the expression on the 
right side of Eq. (3.181) is positive and thus t](s,9) —> +00 as s — » + 00 . The above observations coupled 
with Eq. (3.179) imply that, for a given 0 with 0 < |0| < 7 r, (i) there is one and only one positive value 
such that 

r <0 ifM <v*{o) 

r/ {v 2 ,9) < = 0 if \v\ = v*(9) v ^ 0; 0 < |0| < 7r (3.196) 

{ > 0 if M > v*(9) 


and (ii) 

f=l/2 if |0| = 27 t/3 
^ \ > 1/2 if 0 < |0| < 7r and |0| ^ 2ir/3 


(3.197) 


By using (i) Eqs. (3.196) and (3.197), (ii) Eqs. (3.179) and (3.181), (iii) Theorem 5, (iv) Proposition 2, 
(v) the relation <ti(i/, 0)cr 2 (i / , 9)a^{v 1 9) = —1 (see Eq. (3.22)), and (vi) the fact that the only possible triple 
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root of unit magnitude for Eq. (3.17) is —1 and it occurs only for the case Eq. (3.80) (proved in Sec. 3.6), 
one now arrives at Proposition 3: 

Proposition 3. For a given 9 with 0 < \6\ < i r, Eq. (3.17) has: (i) three distinct roots of unit magnitude 
if 0 < \v\ < i y*(9); (ii) one double root of unit magnitude and a simple root of unit magnitude if \v\ = v*{9)\ 
and (iii) at least one root with its magnitude > 1 if \v\ > 

An immediate result of Eq. (3.197) and Proposition 3 is the following corollary: 

Corollary to Proposition 3. For a given 6 with 0 < \6\ < n, Eq. (3.17) has: (i) three distinct roots 
of unit magnitude if (a) 0 < \u\ < 1/2, and also if (b) \v\ = 1/2 and \9\ ^ 27r/3; (ii) one double root of unit 
magnitude and a simple root of unit magnitude if \v\ = 1/2 and \9\ = 27 t/ 3 ; and (iii) at least one root with 
its magnitude > 1 if \v\ > 1/2 and \9\ = 2tt/3. 

With the above preliminaries, a rigorous study of the stability conditions of the a(3) scheme will be 
given in Sec. 3.9. 

3.9. Stability condition of the a(3) scheme 

Because of Eq. (3.1) and a reason presented right before Sec. 3.1, we have the following definition: 

Definition 1. The a(3) scheme is said to be stable with respect to a given v if, for every 9 , — 7r < 9 < n, 
every element of the matrix [G(i>, 9)] m remains bounded as the positive integer m — * +oo. On the other 
hand, the scheme is said to be unstable with respect to a given v if, for any 9 , — 7r < 9 < 7r, at least one 
element of the matrix \G{v, 9)} m becomes unbounded as m — > +oo. 

To study stability of the a(3) scheme, in the following we introduce needed matrix preliminaries. First 
note that an N x N matrix has at least one eigenvector and at most N linearly independent eigenvectors 
[76]. Related to this subject, we have Definition 2 [76]: 

Definition 2. An N x N matrix A is said to be nondefective if it admits N linearly independent 
eigenvectors. On the other hand, A is defective if it admits less than N linearly independent eigenvectors. 

Another definition we need later is Definition 3: 

Definition 3. Let A^, i = 1,2,3,..., N, be the eigenvalues (which may or may not coincide with one 
another) of an N x N matrix A. Then the spectral radius of A is 

p(A) = f max{|A^|} (3.198) 


Next we have the following well-established theorem [76]: 

Theorem 6. For any N x N matrix, (i) each distinct eigenvalue of multiplicity rn is associated with 
at least one eigenvector and at most m linear independent eigenvectors; and (ii) two eigenvectors associated 
with two different eigenvalues are linearly independent. 

By using Definition 2 along with Theorems 4 and 6, we have Theorem 7, i.e. , 

Theorem 7. An N x N matrix A is defective if and only if A has at least one eigenvalue which satisfies 
the following properties: (i) its multiplicity in is greater than one; and (ii) the number of linearly independent 
eigenvectors associated with this eigenvalue is less than m. 

Moreover, with the aid of Theorems 2 and 3, one can easily prove Theorem 8: 

Theorem 8. Let (i) A be an TV x N matrix, (ii) A be the complex conjugate of A, and (iii) B be a 
matrix similar to A (i.e., there exist a nonsingular N x N matrix S such that B = S~ 1 AS). Then A is 
defective (nondefective) <t=> A is defective (nondefective) <t=> B is defective (nondefective). 

An immediate result of Theorem 7 is Lemma 3, i.e., 

Lemma 3. An N x N diagonal matrix is nondefective. 

Next we will prove Theorem 9, i.e., 
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Theorem 9. Let A be a 3 x 3 matrix. Then every element of A m remains bounded as the positive 
integer m —> +oo if and only if 

/ a\ { < 1 if A is nondefective , 0 lnrA 

' ,( ' 4| Ul if A is defective (3199) 

Proof. Let A be nondefective. Then the Jordan canonical form theorem [76] implies that there is a 
nonsingular 3x3 matrix S so that A = SA 0 S -1 where 


/Ai 0 0 

A 0 d = 0 A 2 0 
V 0 0 A 3 


(3.200) 


with Ai, A2, and A3 being the eigenvalues of A. Because A = SA 0 S _1 implies that A m = S(Ao) m S~\ every 
element of A remains bounded as m — » +00 every element of Ao remain bounded as in — > +00. As such, 

for the nondefective case, Theorem 9 now follow from Eq. (3.198) and the fact that (i) 


(A 0 ) m = 


and (ii) for a complex number c 


m = 1, 2, 3 . 


lim | c 71 

m— >■+ OO 


< 1 if |c| < 1 
= +00 if c > 1 


(3.201) 


(3.202) 


Next let A be defective. Then, according to Theorems 4, 6, and 7, it must belong to one of the following 
mutually exclusive cases: (a) it has an eigenvalue A of multiplicity m = 3 and it admits one and only one 
linearly independent eigenvector; (b) it has an eigenvalue A of m = 3 and it admits two and only two linearly 
independent eigenvectors; and (c) it has an eigenvalue Ai of m = 1 and another eigenvalue A2 of m = 2 
such that there is one and only one linearly independent eigenvectors associated with the eigenvalue A2. 
According to the Jordan canonical form theorem, for each case, there exists a nonsingular 3x3 matrix S 
such that A = SAS _1 where (i) for case (a), 


(ii) for case (b), 


and (iii) for case (c), 


/A 1 0 

A = As = f | 0 A 1 
V 0 0 A 


/ A 0 0 

A = A 2 = f | 0 A 1 
V 0 0 A 


/ Ai 0 0 

A = A 3 = f 0 A 2 1 

V 0 0 A 2 

By induction, one can prove easily the following relations: 

/ A m mA™- 1 [m(m — l)/2]A m-2 \ 
(Ai) m = 0 A m mA”- 1 

\00 A m 


m = 1 , 2 , 3 , 


(3.203) 


(3.204) 


(3.205) 


(3.206) 


(A 2 ) m = 


0 

toA™ -1 

A m 


m = 1 , 2 , 3 , 


(3.207) 
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and 

/(Ai) m 0 0 \ 

(A 3 ) m = 0 (A 2 ) m m(A 2 ) m ” 1 m= 1,2,3,... (3.208) 

V 0 0 (A 2 ) m ) 

Because A = SAS~ l implies that A m = S(A) m S~ 1 , every element of A remains bounded as to — > +oo <+■ 
every element of A remain bounded as to — » +oo. As such, for the defective case, Theorem 9 now follows 
from (i) Eqs. (3.198) and (3.203)-(3.208), and (ii) the fact that Eq. (3.202) and the relations 


lim | me™" 1 1 

m — >-+oo 


0 if |c| < 1 
+oo if jcj > 1 


(3.209) 


and 


lim |[m(m — l)/2]c m 2 | 

m — »+oo 


0 if |c| < 1 

+oo if jcj > 1 


(3.210) 


are true for any complex number c. QED. 

An immediate result of Definition 1 and Theorem 9 is Lemma 4: 


Lemma 4 The a(3) scheme is stable with respect to a given v if and only if, for this v and every 9 , 
— 7T < 9 < 7T, 


p{G[y,9)) 


< 1 if G(v , 9) is nondefective 

< 1 if G(v, 6) is defective 


— 7T < 9 < 7T 


(3.211) 


Because oe(v,9), t = 1,2,3, are the eigenvalues of G(v,9), with the aid of Proposition 1(b) (or part 
(iii) of Corollary to Proposition 3) and Eq. (3.198), Lemma 4 implies that the a(3) scheme is unstable if 
\v\ > 1/2. In the following, we will prove Proposition 4, i.e., 

Proposition 4. The a(3) scheme (i) is stable if and only if 


H < 1/2 (3.212) 

(ii) is neutrally stable for any v satisfying Eq. (3.212); and (iii) is linearly unstable (in a sense to be defined) 
if \v\ = 1/2. 

As a preliminary, first we will study defectiveness of G(v,9) over several domains of iy,6). We begin 
with Lemma 5: 

Lemma 5. (i) G(is, 0) is nondefective for any z/; (ii) G(z^, n) is defective if and only if \v\ = and 

(iii) G(0, 9), —tt < 9 < n, is nondefective. 

Proof. Because 

/ 1 0 0 

G(i/,0) = 0 -1 0 

\0 0 1 

Part (i) follows from Lemma 3 and Eq. (3.213) immediately. 

Next, by using Eq. (3.76), one concludes that G(z^, 7r) has an eigenvalue with multiplicity in > 1 if and 
only if either (i) the two eigenvalues given in Eq. (3.77) are equal, i.e., 

v 2 (2v 2 - 1) = 0 (3.214) 


(3.213) 


or (ii) the eigenvalue —1 is equal to one of those given in Eq. (3.77), i.e., 

1 - 2z^ 2 = v / 2i/ 2 (2z/ 2 - 1) or 1 - 2z/ 2 = -sj2v 2 {2v 2 - 1) 
Eq. (3.215) implies that (1 — 2v 2 ) 2 = 2v 2 {2v 2 — 1) which <+< 

2v 2 = 1 


(3.215) 


(3.216) 
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Combining Eqs. (3.76), (3.214) and (3.216), one concludes that: (i) G(v, tt) has an eigenvalue with m > 1 if 
and only if either v = 0 or v = ±l/\/2; (ii) the eigenvalue of G(0, 7r) with to = 2 is 1; and (iii) the eigenvalue 
of G(±l/v / 2, tt) with to = 3 is —1. Also it can be show easily that (iv) 


1 and 0 


(3.217) 


are two linearly independent eigenvectors of 


/ 3 04/3 

G(0,tt) =01 0 

\ -6 0 -3 


(3.218) 


associated with the eigenvalue 1. By using the above results (i), (ii), and (iv), Theorem 7 implies that G(u , tt) 
is nondefective if \v\ ^ 1 / a/ 2 • To complete the proof of part (ii) of Lemma 5, we need only to show that 
G(±l/y / 2, n) is defective. 

To proceed, note that 


/ 3 t2v / 2 8/3 

G(± 1/n/2,tt)= ±2\/2 -1 ±4^/3 

\ -6 ±3-^/2 -5 


(3.219) 


Let x = (xi,X2,X3 Y be an eigenvector of G(±l/v / 2, n) with the eigenvalue — 1, i.e. , (i) x ^ 0 and (ii) 


By using Eq. (3.219), Eq. (3.220) 


G(±l/-\/2, tt)x = —x 


X 2 = 0 and 3a; i + 2a;3 = 0 


Thus any eigenvector of G(±l/v / 2, tt) associated with the eigenvalue —1 must be in the form 


( 2 

c 0 


(3.220) 


(3.221) 


(3.222) 


where c is a complex number ^ 0. In other words, there is one and only one linearly independent eigenvector 
of G(±l/v / 2, tt) associated with the eigenvalue —1. Because to = 3 for this eigenvalue, Theorem 7 implies 
that G(±l/v / 2, 7r) is defective. Thus the proof of part (ii) is completed. 

Next, according to Eq. (3.75), for the matrix 


/ 2 — cos0 isin0 (2/3) (1 — cos0) 
G(0,9) = I isin# —cos 9 (2/3)isin0 

\3(cos0— 1) — 3isin6* 2 cos 6 — 1 


— 7T < 9 < 7T 


(3.223) 


the eigenvalue with to = 2 is 1 while the eigenvalue with to = 1 is —1. On the other hand, 


i( 1 — cos 9) anc j 27(1 — cos 9) 
sin 9 sin 9 


0 < \9\ < TT 


(3.224) 
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are two linearly independent eigenvectors of G(0, 6 ) associated with the eigenvalue 1 if 0 < \ 6 \ < n. By using 
the above results, part (iii) of Lemma 5 now follows from Theorem 7 and parts (i) and (ii) of Lemma 5. 
QED 

Next, defectiveness of the matrix G(v,9) when \u\ = 1/2 and |0| = 27r/3 will be established in Lemma 
6, i.e., 

Lemma 6. When \v\ = 1/2 and |0| = 2-7 t/3, G(y, 6 ) is defective, and it has an eigenvalue with 
multiplicity m = 2 and another eigenvalue with m = 1. 

Proof. Consider the case v = 1/2 and 9 = 27t/ 3. Eq. (3.62) implies that 


cr 0 (27r/3) =ct_(2tt/ 3) = -(l-n/3)/2 and cr+(27r/3) = (1 - n/3)/2 (3.225) 


i.e., the eigenvalue of G(l/2, 27 t/ 3) with m = 2 is — (1 — iy/2>)/2 while the eigenvalue with m = 1 is (1 — is/Z)/2. 
Moreover, let x = (xi,X 2 ,X 3 y be an eigenvector of 


G( 1/2, 2tt/3) 


/ (10 — i\/3)/4 (—12 + t5 v / 3)/8 

(3 + n/3)/2 (-l + iVS)/4 

\ -9/2 (9-*6-v/3)/4 


(12 — i3 v / 3)/8 \ 
(3 + ij 3)/4 
(-11 + *2\/3)/4 / 


with the eigenvalue —(1 — iv / 3)/2, i.e., (i) x ^ 0 and (ii) 


G(l/2, 2tt/3)x = 


(1 - iV3)/2 


(3.226) 


(3.227) 


By using Eq. (3.226), Eq. (3.227) <t=> 

12-i3\/3 —12 + *5-^/3 \ / x 

3 + i-s/3 1 - i-v/3 ( 2:Cl + 2:3 ) = 0 (3.228) 

3 -3 + Z2V3/V 3:2 ) 


Because the 3x2 coefficient matrix in Eq. (3.228) is formed by two linearly independent 3x1 column 
matrices, Eq. (3.228) <t=> 

2x i + X 3 = 0 and X 2 = 0 (3.229) 

Thus any eigenvector of G(l/2, 27r/3) associated with the eigenvalue —(1 — i\J 3)/2 must be in the form 


c 



(3.230) 


where c is a complex number ^ 0. In other words, there is one and only one linearly independent eigenvector 
of G(1/2,27 t/ 3) associated with the eigenvalue —(1 — i\J 3)/2. Because m = 2 for this eigenvalue, Theorem 
7 implies that G(l/2,27r/3) is defective. 

For each of the matrices G(l/2, — 27 t/3), G(— 1/2, 27t/ 3), and G(— 1/2, — 27r/3), its defectiveness can be 
proved by using (i) defectiveness of G(l/2, 27 t/ 3), (ii) Eqs. (3.4) and (3.5), and (iii) Theorem 8. Also the 
fact that each has an eigenvalue with m = 2 and another eigenvalue with m = 1 can be proved by using 
Eqs. (3.62) and (3.63) directly, or by using Eq. (3.14) along with the proved similar property of G(l/2, 27r/3). 
QED. 

Next, an immediate result of Theorem 7 and part (i) of Corollary to Proposition 3 is Lemma 7: 
Lemma 7. For a given 9 with 0 < |0| < 7r, G(v,9) is nondefective if (a) 0 < \v\ < 1/2, and also if (b) 
\v\ = 1/2 and \9\ ^ 27r/3. 

Moreover, by combining Lemmas 5 and 7, we have Lemma 8: 
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Lemma 8. G(y,9) is nondefective if (a) \u\ < 1/2 and — tt < 9 < tt, and also if (b) \v\ = 1/2, 
— 7 r < 9 < tt, and \9\ ^ 2tt/3. 

To prove Proposition 4, note that it has been shown earlier that the a(3) scheme is unstable if \v\ > 1/2. 
On the other hand, by using Proposition 1(a) and Definition 3, one has 

\p(G(u,9)) | = \ai(v,0)\ = |<7 2 (z',6>)| = \a 3 (v,9)\ = 1, -tt < 9 < tt, if \v\ < 1/2 (3.231) 

In turn, with the aid of Eq. (3.231) and part (a) of Lemma 8, Lemma 4 implies that the a(3) scheme is 
neutrally stable if \u\ < 1/2. Thus, to complete the proof, one needs only to prove that the a(3) scheme is 
linearly unstable if \v\ = 1/2. 

By using Proposition 1(a), part (b) of Lemma 8, and Theorem 9, one concludes that every element of 
[G(u,9)] m remains bounded as to — > +oo for any with \u\ = 1/2, —tt < 9 < tt, and \9\ ^ 2tt/3. Thus, 
Definition 1 implies that the a(3) scheme is unstable at v = ±1/2 if and only if at least one element of 
[G(±l/2, 0 )] m becomes unbounded as m — > ±00 when 9 = 2tt/3 or 9 = —2tt/3. 

Consider any ( v,9 ) with \u\ =1/2 and \9\ = 2tt/3. Then, by using Eq. (3.231) and Lemma 6 , Theorem 
9 implies that at least one element of [G(is, 9)] m becomes unbounded as to — > ± 00 , i.e., we have proved that 
the a(3) scheme indeed is unstable if \v\ = 1/2. Moreover, according to Lemma 6, G(y,9) is defective and 
has an eigenvalue (denoted by ct\{v,9)) with multiplicity = 1 and another eigenvalue (denoted by ct 2 (v, 9 )) 
with multiplicity = 2. As such the Jordan canonical form theorem implies that there exists a nonsingular 
3x3 matrix S such that G{v, 9) = SA 3 S~^ where A 3 is defined in Eq. (3.205) with 

Ai = tri(^,0) and A 2 = 02 {it, 9) (\u\ = 1/2* \ 6 \ = 2tt/Z) (3.232) 

Because G(v, 9) = SA. 3 S - 1 implies that [G(iy 9)] m = S(A 3 ) rn S~ 1 , the behavior of the elements of [G{v, 9)] m 
as to — > ±00 is completely determined by that of (A 3 ) m as to — > ± 00 . Moreover, because |Ai| = | A 2 I = 1 
(which follows from Eqs. (3.231) and (3.232)), Eq. (3.208) implies that the only element of (A 3 ) m that will 
become unbounded as m — > ±00 is the element m( \ 2) m ~ 1 and that the magnitude of this element grows 
only linearly with to. As a result of the above considerations, one concludes that any element of [G(y, 0)] m 
at most can only grow linearly with to for any case with \v\ = 1/2 and \9\ = 2n/3. Because of this reason 
and the fact that the time evolution of the round-off errors originally introduced during any marching step 
is also governed by the a(3) scheme, the round-off errors originally introduced during a single marching step 
at most can only grow linearly with the marching-step number if \v\ = 1/2. It is in this sense that the a(3) 
scheme is said to be linearly unstable when \v\ = 1/2. QED. 

Note that the total round-off errors observed at the start of any marching step is the sum of the 
“offsprings” of the round-off errors originally introduced during all the previous marching steps. Because 
of the intrinsic random nature of round-off-error generation and the accompanied effect of (in-phase and 
out-of-phase) interferences, evaluating the sum referred to above is much more complex than evaluating the 
offspring of the round-off errors introduced during a single marching step. Nevertheless, because the growth 
rate of the magnitude of the term m(\ 2) rn ~ 1 for the case | A 2 1 = 1 is much lower than the exponential growth 
rate of a term associated with a case with p(G(y,9)) > 1, one still can infer that the instability of the a(3) 
scheme when \n\ = 1/2 is much milder than that for a case with \v\ > 1/2. As will be shown in Sec. 4, this 
prediction is borne out by numerical results. 
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4. Numerical results 


To assess the accuracy of the a(3) scheme, consider the model problem with the PDE 


du 

at 



(4.1) 


and the exact solution 


u = u e (x,t) = sin [27r(x 





i2-K(x—t) 


e -i2ir{x-t) 


(4.2) 


We have 

a = L = T = 1 (4.3) 

where L = wavelength and T = period. Moreover, u e (x, t ) is a linear combination of two plane wave solutions 

e k + (x-i) and e k.(x-t) with 

k± = ±2 tt (4.4) 


Let (i) 


, ., def du e (x,t) 

u xe (x,t) = — 


and u xxe (x,t) d = ^ U ' 


(x,t) 


dx 2 


(4.5) 


and (ii) the spatial domain of unit length be divided into K uniform intervals. Then, with the aid of Eq. (4.3), 
one has 

Ax = 1 /K, At = is ax, and t = n&t (4-6) 


where n = number of time steps, and t = total marching time. The computer code solving the model 
problem for various values of K and n using the a(3) scheme is listed in Appendix C while that using the 
dual a scheme [71] is listed in Appendix D. Because the dual a scheme (which are defined over the set ft) 
is formed by two completely decoupled a schemes (which are defined over the sets fli and respectively), 
the accuracy and stability conditions of the dual a scheme are identical to those of the a scheme. 

Based on the von Neumann analysis, it was shown in Sec. 3 (Proposition 4) that the a(3) scheme (i) is 
stable if and only if \v\ < 1/2; (ii) is neutrally stable if \v\ < 1/2; and (iii) is linearly unstable if \v\ = 1/2. On 
the other hand, by using the amplification matrix given in Eq. (3.51) of [71] and a line of arguments similar 
to that presented in Sec. 3, one can show that the a scheme (and the dual a scheme) (i) is stable if and only 
if \v\ < 1; (ii) is neutrally stable if \v\ < 1; and (iii) is linearly instable if \v\ = 1. These stability conditions 
have been verified numerically using the codes implementing the a(3) scheme and the dual a scheme, which 
are listed in Appendices C and D, respectively. 

In Tables 1-4, the numerical errors of several computations using the a(3) scheme and the dual a scheme 
are presented in terms of the following error norms for the non-normalized independent mesh variables: 


E(K, n, v) 


def 


N 


l 

K 


K - 1 


j—0 


U e (Xj,t n )] 2 


and 


E x (K,n,v) = f 


\ 


1 K ~ X 

J7 u xe{Xj,t n )\ 2 

3=0 


E„(K, n,v) = 


\ 


K - 1 


~ Uxxe(Xj,t n )] 2 

3=0 


(4.7) 


(4.8) 


(4.9) 
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Because, at each mesh point ( j,n ), the only non-normalized independent mesh variables associated with 
the dual a scheme are w” and (u x )”, obviously the error norm E xx (K,n,u) is not applicable to the dual a 
scheme. 

The numerical errors of several simulations with v = 0.1 and t = 9.876 are given in Table 1. For the 
dual a scheme, as the values of K and n become larger, the values of E and E x are both reduced by a factor 
of about 4 as both K and n double their values, i.e. , the scheme is 2nd order in accuracy for both it” and 
( u x )”. On the other hand, for the a(3) scheme, the values of E, E x , and E xx are reduced by the factors 
16, 16, and 4, respectively. Thus the a(3) scheme is 4th order in accuracy for both it” and (w x )” while only 
2nd order in accuracy for (it xx )”. From the results shown, one can see that the a(3) scheme is much more 
accurate than the dual a scheme. As an example, for the case with K = 25 and n = 2469, the value of E 
for the dual a scheme is larger than that for the a(3) scheme by a factor of 3450! Because the a scheme is 
only 2nd order in accuracy for it”, it is estimated that the accuracy of it” achieved by the a(3) scheme with 
K = 25 and n = 2469 is identical to that achieved by the dual a scheme with K = 25 x ^3450 ss 1468 and 
n = 2469 x V3450 « 145029. 

Here the reader is reminded that, for a reason given in Sec. 2.5, the conclusion reached above about 
the orders of accuracy of the a(3) scheme does not contradict that reached in Sec. 2.5 about the orders of 
truncation error of the a(3) scheme. 

In Table 2, the cases considered have v = 0.1 and t = 10.00 = 10T. For these cases where t is an integer 
multiple of the period T, it is seen that the a(3) scheme is 4th order in accuracy for it”, (its,)”, and (it xx )”. 

In Table 3, the cases considered have v = 0.5 and t = 49.38, For these cases where the value of v is 
right at the stability boundary of the a(3) scheme, aside from round-off errors, the numerical values of (u x )” 
generated using the a(3) scheme are all identical to their exact solution values, respectively, if n are even 
integers. 

In Table 4, the cases considered have v = 0.5 and t = 50.00 = 50T, i.e., the value of v is right at the 
stability boundary of the a(3) scheme and t, is an integer multiple of T . It is seen that the numerical values 
of it" , (tt x )”, and (it xx )” generated using the a(3) scheme, aside from round-off errors, are all identical to 
their exact solution values, respectively. Note that: (i) n and At are chosen according to Eq. (3.67) and n is 
even for each of these cases, and (ii) the exact solution are a linear combination of two plane wave solutions 
with \6\ = |fc±Ax| = | ± 2tt/K\ < 27r/3, K = 25,50,100,200 (see Eq. (4.4)), i.e., 0 observes the condition 
Eq. (3.64). Thus the numerical results of the a(3) scheme shown in Table 4 confirm an accuracy prediction 
made in Sec. 3.4. 

Moreover, the round-off errors associated with the a(3) scheme shown in Table 4 also confirm a prediction 
made at the end of Sec. 3.9, i.e., the a(3) scheme is only mildly unstable when \u\ = 1/2. 

According to the von Neumann analysis, (i) the dual a scheme has no dissipative errors (i.e., the 
magnitudes of all its amplification factors = 1 for all 9 in the domain — 7r < 9 < tt) if \v\ < 1; and (ii) the 
a(3) scheme has no dissipative errors if \v\ < 1/2. Thus, for a simulation with periodic boundary condition, 
aside from round-off errors, the numerical errors generated using the dual a scheme are contributed solely 
by the phase (dispersive) errors if \v\ < 1. On the other hand, those generated by using the a(3) scheme are 
contributed solely by the phase errors if \v\ < 1/2. Thus the relative accuracy of the a scheme and the a(3) 
scheme can also be evaluated by comparing their phase errors. 

For both the a(3) scheme and the dual a scheme, the phase error of the principal amplification factor 
associated with any (y 1 9) can be measured by 

E P (u, 9) = f 1 - (M^ 0) ^ 0) (4.10) 

<Pe{V,9) 

Here: (i) (j> e (v,9) is the phase angle of the exact amplification factor, i.e., 

(j>e{v, 9) = -vQ (4-11) 

and (ii) <f>(u,9) is the phase angle of the principal amplification factor. Note that, by using Eq. (3.14) and 
similar relations for the dual a scheme, one concludes that 

9) = -9) = 9) = -</(-ji, -9) (4.12) 
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An immediate result of Eqs. (4.10)-(4.12) is 

E p (v, 9) = E p (-u , 9) = E p (u, -9) = E p (—v, -9) (4.13) 

Thus only the values of E p (u,9) with nonnegative v and nonnegative 9 need to be evaluated numerically. 

In the code listed in Appendix E, for the dual a scheme, </>(r/, 9) is evaluated using the exact formula: 

v 2 < 1 — tv < 9 < n (the dual a scheme) (4-14) 

(see Eq. (3.31) in [71]). On the other hand, for the a( 3) scheme, <j)(u, 9) is evaluated in the same code using 
the Newton’s iterative procedure Eq. (3.113) and the assumption 

4>° = <j) e (y,9) (4-15) 

After k iterations, <p k is taken as the converged value of <j>(v, 9) if 

\{4 >k / 4 >k ~ 1 ) — 1 1 < e (4.16) 

where e > 0 is a very small preset value. Note that the iterative procedure generally converges rapidly. 
Specifically, Eq. (4.16) with e = 10 -14 is satisfied after at most 5 iterations for all (y,9), \v\ <1/2 and 
— 7 r < 9 < n. Moreover, as expected, convergence is reached after one iteration for all 9 if |^| = 1/2. 

The numerical values of E p {v,9) for the cases v = 0.001,0.01,0.1,0.5 are plotted against 9 (denoted by 
Z) in Fig. 3 for both the dual a scheme and the a(3) scheme. The values of Ep(^, 9) for the dual a scheme are 
calibrated using the left-ordinate scale while those for the a(3) scheme are calibrated using the right-ordinate 
scale. It can be seen that the values of E p (v, 9) for the a(3) scheme are uniformly much smaller than those 
for the dual a scheme. In fact, the numerical results indicate that, for the a(3) scheme, (i) cj)(y,9) = O [(0) 4 ] 
if |^| < 1/2; and (ii) aside from round-off errors, 0(^,0) = 0 for all 9 , if |^| = 1/2 (which is expected from 
Eqs. (3.62) and (3.63) — see a discussion given following Eq. (3.64)). On the other hand, for the dual a 
scheme, c/)(v,9) = O [(0) 2 ] if |^| < 1. 


</>(^, 9) = tan 


— j/sin 9 
— v 2 sin 2 9 
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5. Conclusions and discussions 


A thorough and rigorous discussion of the a(3) scheme, a new high order neutrally stable CESE solver 
of Eq. (1.1) has been presented. As in the case of other similar CESE neutrally stable solvers [1,5,11,72], 
the a(3) scheme enforces conservation laws locally and globally, and it has the basic, forward marching, 
and backward marching forms. These forms are equivalent and satisfy the PT invariant property defined in 
Sec. 2. 

Based on the concept of PT invariance, the algebraic relations Eqs. (2.114)-(2.118) are derived in Sec. 2. 
As it turns out, in the von Neumann analysis presented in Sec. 3, these relations can be used to prove that 
the a(3) scheme is neutrally stable when it is stable. Another set of algebraic relations Eq. (2.119) which 
results from other invariant property are also discussed in Sec. 2. 

In addition to establishing the neutral stability of a(3) scheme, in Sec. 3, it is also proved rigorously 
that all three amplification factors of the a(3) scheme are of unit magnitude for all phase angles 9 if and 
only if \v\ < 1/2 (Proposition 1). Moreover, it is proved that the a(3) scheme is (i) stable if and only if 
\v\ < 1/2; and (ii) linear unstable if \v\ = 1/2 (Proposition 4). These theoretical results have been confirmed 
by numerical experiments. 

It is shown in Sec. 4 that the a(3) scheme generally is (i) 4th-order accurate for the mesh variables 
u" and (u x )j”, and (ii) 2nd-order accurate for (u^)". However, in some exceptional cases, the scheme can 
achieve perfect accuracy aside from round-off errors. Moreover, it is shown that the phase errors of the 
principal amplification factor of the a(3) scheme are 0(9 4 ) if \v\ < 1/2, a sharp reduction from those of the 
dual a scheme which are 0(9 2 ) if \v\ < 1. 
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Appendix A. Proof for Theorems 1 and 2 


First we prove Theorem 1. According to Eq. (8.20) on p.265 of [75], the fact that Xe, £ = 1,2 ,N 
are the eigenvalues of the N x N matrix A ■<=> 

det(A - XI) = (Ai - A)(A 2 - A) • • • (X N - A) (A. 1) 

where A is any complex variable and / is the N x N identity matrix. Let A = 0. Eq. (A.l) implies that 


det (A) = AiA 2 • • • A tv 


(A2) 


By definition, A is nonsingular •<=>■ det(A) 0. Thus part (i) follows from Eq. (A. 2). 
According to Eq. (A.l), to prove part (ii) we need only to show that 


det(A- 1 -A/) = (--A)(--A)...(--A 


1 


A 2 


1 


A N 


(A3) 


for any complex variable A. Because det (BC) = det(B) det(C) for any two N x N matrices B and C, we 
have det(A)det (A -1 ) = det (AA _1 ) = det(/) = 1, i.e., det(A _1 ) = l/det(A). Thus Eq. (A. 2) implies that 


det (A 1 ) = 


1 


AiA 2 • • • Ajv 


By comparing Eqs. (A. 3) and (A. 4), one concludes that Eq. (A. 3) is valid if A = 0. 
Let A 0. Then 

1 


Thus 


A -1 - XI = -A I A- 1 ( A - A/ 

A 


clet (A 1 — A/) = det(— XI) det (A 1 ) det ^A — —I^j 


(A.4) 


(A.5) 


(A.6) 


With the aid of (i) Eqs. (A.l) and (A.4) and (ii) the fact that det(— XI) = (— A) , Eq. (A.6) implies that 


det (A 1 — A I) = 


(—A) 


JV 


A1A2 • • • Xn 


Ai - t ) A 2 - - 


Xn ~x ] = 


-A 


Ai 


A 


-A 


A 2 


A 


-A 


Xn 


Ar — - — — A 2 — T • • • - Ajv— T = _ r — A { — A • • • { A 


Ai 


A 2 


Xn 


(A. 7) 


i.e., Eq. (A. 3) is also valid if A 0. Thus part (ii) of Theorem 1 has been proved. QED. 
According to Eq. (A.l), to prove Theorem 2, we need only to show that 

det (A - XI) = (AT - A)(A 2 - A) • • • (Xn - A) 


Because clet(Af) = det (AI), by using Eq. (A.l), we have 

clet(A — XI) = det^A — A/j = det(A — XI) 


(A. 8) 


(A.9) 


= (Ai - A)(A 2 - A) • • • (Ajv - A) = (A, - A)(A 2 - A) • • • (X N - A) 


i.e., Eq. (A. 8) has been proved. QED. 
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Appendix B. Proof for Eqs. (3.138), (3.139), (3.155), (3.169), and (3.170) 

Proof for Eq. (3.138). Assuming Eq. (3.134) and using elementary trigonometry, we have 


cos 0i -(- cos 02 — cos(0i -(- 0 2 ) — 1 _ 2cos ++ 2 )cos + 1 + 2 ) — 2cos 2 ++ 2 ) 

sin cj>i + sin 0 2 + sin(0i + 0 2 ) ~~ 2 sin(^±^) cos(^^) + 2 sin(^±^) cos(^±^) 


cos( 


sin( 


0 1+02 
2 


01+02 

2 


) cos( 
) cos( 


^=^) -cos( 
+ cos( 


01+02 \ 
2 ' 


01+02 \ 
2 1 


= cot 


01 + <t> 2 


2 sin(0i /2) sin(0 2 /2) 
2cos(0i/2) cos(0 2 /2) 


= cot 


^2 ^ 


tan(0i/2)tan(0 2 /2) = tan(0i/2) tan(0 2 /2) tan ( ±— — 


7T 01+02 


Eq. (3.138) follows from Eq. (B.l) and Eq. (3.130). QED. 

Proof for Eq. (3.139). Assuming Eq. (3.134) and using elementary trigonometry, we have 


(B.l) 


cos 0i + cos 0 2 - cos(0i + 0 2 ) + 3 = 2cos(&±&) cos(^U^) + 2cos 2 (^+ + 4sin 2 (^±^) 
sin + sin 02 + sin(0i + 0 2 ) 2sin(+t^) cos C^ 1 "^ 2 ) + 2sin ++ <fe ) cos +)+ ) 


cos( 01 + 2 ) 

[cOS+ 1 + 2 ) + COS+ 1 + 2 )l 

+ 2sin 2 + + 2 ) 

( + + 02 ^ 

, sin(^) 

sin(^+^) 

[cos( 01 " 02 ) + cos++ 02 )] 

U | 

l 2 J 

cos(0i/2)cos(0 2 /2) 

(B.2) 


/ 7T 0i + 0 2 \ sin(0i/2)cos(0 2 /2) +cos(0i/2)sin(0 2 /2) 

an V 2 2 ) cos(0i/2)cos(0 2 /2) 


^ + tan(0!/2) + tan(0 2 /2) 


Eq. (3.139) follows from Eq. (B.2) and Eq. (3.130). QED 

Proof for Eq. (3.155). Using elementary trigonometry, we have 


i • i • / , \ r, • (Pi + 02^ (Pi ~ 02^ 0 . /+1 + 02^ /+1 + +2n, 

sinipi + sin <p 2 + sm((p! + </j 2 ) = 2sm( ) cos( ) + 2sm( ) cos( ) 

0 . (Pl+P2,\ (Vl-Vls . (Pl + P2P\ _ . . ( Pl+P2, , , /0 S 

= 2 sm( ) j^cos( ) + cos( )J = 4sm( ) cos(<^ 1 /2) cos(<0> 2 /2) 

Next, by using Eq. (3.153), we have 


(B. 3) 


. (P1+P2, . /,7T +3\ , /0 x 

sm( ) = sin \ ±- - —j = ±cos(0 3 /2) 


(BA) 


Eq. (3.155) is a direct result of Eqs. (B.3) and (B.4). QED. 

Proof for Eq. (3.169). Using elementary trigonometry, we have (i) 


tan(0/4) — cot(0/4) = tan(0/4) — 


1 


tan(0/4) 


2[1 — tan 2 (0/4)] 
2 tan(0/4) 


tan(0/2) 


0 < \6\ < 7T (B. 5) 
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and (ii) 


0 < \e\ < 7T (B.6) 


0 < \0\ < 7T {B. 7) 
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2 //w „, 2[1 + tan 2 (0/2)l 2sec 2 ((9/2) 

- tan(0/2) = tan(0/2) = tan(0/2) 

tan(0/2) tan(0/2) tan(0/2) 

2 4 

= tan(0/2) 777 7— . 7- 7— = tan(0/2) ; — - 

v ' cos(0/2) sin (0/2) v ' ’ sin0 

Eq. (3.169) is a result of Eqs. (B.5) and (B.6). QED. 

Proof for Eq. (3.170). Using elementary trigonometry, we have 

tan(0/2)[cot(0/4) - tan(0/4)] = tan(6>/2) 1 = tan(0/2) = 2, 


i.e., Eq. (3.170) has been proved. QED. 



c Appendix C. Code "a3.for'' 

c 

implicit real*8 (a-h,o-z) 
parameter (nxd=2000) 

dimension u(nxd), un(nxd), ux(nxd), uxn(nxd), uxx(nxd), 

* uxxn(nxd) 
c 

c Code "a3.for“. This code implements a neutrally stable solver 

c for a pure advection equation with a sine-wave initial data and 

c periodic boundary conditions. The sine wave is of unit wavelength, 

c 

c Theoretically, the solver is designed to have at least 

c third-order accuracy. However, it has been shown numerically 

c that the scheme is of 4th-order accuracy, 

c 

c Let a, dx, and dt denote the advection speed, the spatial mesh 
c interval and the time step size, respectively. Let (i) the 

c Courant number cn = a*dt/dx; and (ii) theta denote the phase 

c angle of a Fourier mode. Then, according to Proposition 1(a), 

c for any pair of cn and theta with |cn| . le. 1/2 and -pi .It. 

c theta .le. pi, the three amplification factors of the a3 scheme 

c are all of unit magnitude, 
c 

c It can be shown analytically the scheme is (i) stable if 

c jcn| .It. 1/2; (ii) linear unstable if |cn| = 1/2 ; and 

c (iii) unstable if |cn| .gt. 1/2. 

c 

c There are three independent mesh variables (the analogues of the 

c dependent variable and its first and second spatial derivatives) 

c and three eauations per mesh points. The three equations are 

c obtained by imposing (i) two conservation conditions over CE- and 

c CE+ and (ii) a STI invariant condition that insures the solver has 

c third-order truncation errors . 

c 

c ux, uxx, uxn, and uxxn represent the current and updated numerical 

c analogues of normalized spatial derivatives during time marching, 
c However, the output contains the non-normalized values, 

c 

c it = no. of time steps, 
c k = mesh intervals per unit length, 

c cn = Courant number . 

c a = advection speed. 

c iop = output selector . Only the global L2 error norms eru2 , erux2 , 

c and eruxx2 (which correspond to u, ux, and uxx, respectively) 

c along with the problem defining parameters will be included in 

c the output if iop = 0. Otherwise, all local solution and error 

c values will also be included, 

c 

c The computational domain is 0 .le. x .le. 1. In the output, (i) 
c ue, uxe and uxxe are local exact solution values, (ii) u, ux, 

c and uxx are local numerical solution values, (iii) eru = u - ue, 

c erux = ux - uxe, and eruxx = uxx - uxxe, and (iv) eru2, erux2 , 

c and eruxx2 (which are given at the end) are the global L2 error 

c norms for u, ux, and uxx, respectively, 

c 

it = 6000 
k = 120 
cn = 0 . 5d0 
a = l.dO 
iop = 0 
c 

dx = l.dO/dfloat (k) 
dt = cn*dx/a 
dxl = dx/2 .d0 
dx2 = dx**2/4.d0 
t = dt*dfloat (it) 
kl = k + 1 
k2 = k + 2 

pi = 3.1415926535897932d0 

tp = 2.d0*pi 

tps = tp**2 

tpat = tp*a*t 

pdx = pi*dx 

pdxs = pdx* *2 

cns = cn**2 

cnp = 1 . dO + cn 
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cnm = l.dO - cn 
al = (1 .dO + 2 ,dO*cns) /3 . dO 

a2 = 2 . dO * ( cnp + cns ) / 3 . dO 
a3 = 2. dO* (cnm + cns)/3.d0 
bp = cnp/2.d0 
bm = cnm/ 2 . dO 
cz = 2.d0*cn 
cm = 0 . 5d0 - cn 
cp = 0 . 5d0 + cn 
c 

open (unit=8 , file= ' a3 . opt ' } 

write (8,10) 

write (8,15) 

write (8,20) it,k,dt 

write (8,25) a,cn 

write (8,30) t 

do 100 j = 1 , k2 

tpx = tp*d£loat ( j-1) *dx 

u ( j ) = dsin(tpx) 

ux(j) = pdx*dcos (tpx) 

uxx(j) = -pdxs*dsin (tpx) 

100 continue 

do 400 i = l,it 
do 200 j = 2 , kl 

s = u(j) - cn*ux(j) + al*uxx(j) 
sp = u(j+l) - cnp*ux(j+l) + a2*uxx(j+l) 
sm = u(j-l) + cnm*ux(j-l) + a3*uxx(j-l) 
un(j) = 2.d0*s - bp*sp - bm*sm 
uxn(j) = cz*s + cm*sp - cp*sm 
uxxn(j) = 1.5d0*(sp + sm) - 3.d0*s 
200 continue 

do 300 j = 2 , kl 
u ( j ) = un ( j ) 
ux ( j ) = uxn ( j ) 
uxx { j ) = uxxn ( j ) 

300 continue 

u(k2) = u(2) 
ux(k2) = ux(2) 
uxx(k2) = uxx (2) 
u ( 1 ) = u { kl ) 
ux(l) = ux (kl) 
uxx(l) = uxx(kl) 

400 continue 

eru2 = O.dO 
erux2 = O.dO 
eruxx2 = O.dO 
do 500 j = l,k 
x = dfloat ( j-1) *dx 
tpx = tp*x 
ue = dsin(tpx-tpat) 
uxe = tp*dcos (tpx-tpat) 
uxxe = ~tps*ue 
ux ( j ) = ux ( j ) /dxl 
uxx(j) = uxx(j)/dx2 
eru = u ( j ) - ue 
erux = ux ( j ) - uxe 
eruxx = uxx ( j ) - uxxe 
eru2 = eru2 + eru* *2 
erux2 = erux2 + erux* *2 
eruxx2 = eruxx2 + eruxx* *2 
if (iop.eq.0) goto 500 
write (8,35) x 
write (8,40) ue,u(j),eru 
write (8,45) uxe, ux (j ), erux 
write (8,50) uxxe , uxx ( j ), eruxx 
500 continue 

eru2 = dsqrt (eru2/dfloat (k) ) 
erux2 = dsqrt (erux2/df loat (k) ) 
eruxx2 = dsqrt (eruxx2 /dfloat (k) ) 
write (8,15) 

write (8,55) eru2 , erux2 , eruxx2 
close (unit=8) 

10 format (' OUTPUT FOR CODE A3') 

^ ^ £ onticit { 1 *★******■*** + ****★***★★* + ****** 

20 format (’ it =',i8,' k =',i8,' dt =',gl4 

25 format {' a =’,gl4.7,' CFL =',gl4.7) 
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30 

format 

< ' t =' 

, g!4 .7) 


35 

format 

C x =' 

,g!4. 7,' ******************************************') 

40 

format 

( ' ue = 

' , gl4 . 7 ,xx. 

' u = 1 , gl4 . 7 , xx, 1 eru =',gl4.7) 

45 

format 

( ' uxe 

= ' ,gl4 .7 ,x, 

1 ux =',gl4.7,x,' erux =',gl4.7) 

50 

format 

( ' uxxe 

=- ,gl4.7, 1 

uxx = ' , gl4 . 7 , ' eruxx = ' , gl4 . 7 ) 

55 

format 

stop 

end 

( 1 eru2 

=• ,gi4.7, ■ 

erux2 =',gl4.7,' eruxx2 =',g!4.7) 
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c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


100 


200 


300 


Appendix D. Code "a2.for u 

implicit real*8 (a-h, o~z) 
parameter (nxd=2000 } 

dimension u(nxd) , un(nxd), ux(nxd), uxn(nxd) 

Code xx a2.for’’. This code implememts the dual ''a 1 ' scheme 
with a sine-wave initial data and periodic boundary conditions . 
The sine-wave is of unit wavelength. It can be shown that 
both exact and numerical solutions must be spatially periodic 
with unit wavelength too. 

ux and uxn represent the numerical analogues of normalized 
spatial derivatives during time marching. However , the output 
contains the non-normal i zed values . 

***************************************************************** 
input : 

it - no. of time-marching steps, 
k = no. of spatial intervals per unit length, 
cn - Courant number, 
a = the advection speed. 

iop - output selector. Only the global L2 error norms eru2 and 

and erux2 (which correspond to u, and ux, respectively) along 
with the problem defining parameters will be included in 
the output if iop = 0. Otherwise, all local exact solution 
and numerical solution values {ue, uxe, u, and ux) along with 
the error values (eru and erux) will also be included. 
+****★***★**★★★★***★★***★******+********★**********■*******★★****★ 

The computational domain is 0 .le. x .le. 1, In the output, (i) 
ue and uxe are local exact solution values, { ii ) u and ux are 
local numerical solution values, (iii) eru = u - ue and 
erux = ux - uxe, and (iv) eru2 and erux2 (which are given at the 
end) are the global L2 error norms for u and ux, respectively. 

it = 20000 
k = 200 
cn = 0 . 5d0 
a = l.dO 
iop = 0 

pi = 3. 1415926535897932 dO 

dx = 1 .d0/dfloat (k) 

dt = cn*dx/a 

hdx = dx/2 .d0 

t = dt*df loat (it) 

kl = k + 1 

k2 = k + 2 

tp = pi*2 . dO 

tpat = tp*a*t 

pdx = pi*dx 

cns = (l.dO - cn**2)/2.d0 
cnp = (l.dO + cn)/2.d0 
cnm = (l.dO - cn)/2.d0 

open (unit=8 , file= ' a2 . opt 1 ) 

write (8,10) 

write (8,15) 

write (8,20) it,k,dt 

write (8,25) a,cn 

write (8,30) t 

do 100 j = 1, k2 

tpx = tp*dfloat ( j-1) *dx 

u(j) = dsin(tpx) 

ux(j) = pdx*dcos (tpx) 

continue 

do 400 i = l,it 

do 200 j = 2 , kl 

un(j) = cnm*u(j+l) + cnp*u(j-l) + cns*(ux(j-l) - ux(j+l)) 

uxn(j) = (u(j+l) - u(j-l))/2.d0 - cnp*ux(j+l) - cnm*ux(j-l) 

continue 

do 300 j = 2 , kl 

u(j) = un(j) 

ux ( j ) = uxn ( j ) 

cont inue 
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400 


u(k2) = u (2 ) 
ux(k2) = ux ( 2 ) 
u(l) = u (kl) 
ux(l) = ux(kl) 
continue 
eru2 = 0 . dO 
erux2 = 0 . dO 
do 500 j = 1 , kl 
x = df loat ( j -1) *dx 
tpx = tp*x 
ue = dsin(tpx-tpat) 
uxe = tp*dcos (tpx-tpat) 
ux(j ) = ux( j ) /hdx 
eru = u(j) - ue 
erux - ux ( j ) - uxe 
eru2 = eru2 + eru* *2 
erux2 = erux2 + erux* *2 
if (iop.eq.O) goto 500 
write (8,15) 
write (8,35) x 
write (8,40) ue,u (j ) , eru 
write (8,45) uxe,ux( j ) , erux 
500 continue 

eru2 = dsqrt (eru2/dfloat (k) ) 
erux2 = dsqrt (erux2/dfloat (k) ) 
write (8,15) 
write (8,50) eru2 , erux2 
close (unit=8) 


10 

format 

( 1 

OUTPUT FOR CODE 

A2 1 ) 


15 

format 

( 1 

************************************* 

***************** 

20 

format 

( ' 

it = 

, i8 , 1 k = 1 

i8 , 1 dt =• ,gl4,7) 


25 

format 

( ’ 

a = 1 

gl4 . 7 , 1 CFL =' , gl4 .7) 


30 

format 

(' 

t - 1 

gl4 .7) 



35 

format 

C 

x = ’ 

gl4.7) 



40 

format 

( ’ 

ue - 

, gl4 .7, xx. 

u = ' , gl4 . 7 ,xx, ' eru 

=' ,gi4.7) 

45 

format 

(’ 

uxe 

, gl4 . 7 , x, 

ux = 1 , gl4 . 7 , x, ' erux 

=' ,gl4.7) 

50 

format 

(’ 

eru2 

=' ,gl4.7, ’ 

erux2 = ' , gl4 . 7 ) 



stop 

end 
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Appendix E . Code " f a3 . for 


c 
c 

implicit real*8 {a-h, o-z) 
c 

c Code ' ' fa3 . for' ' 
c 

c Let a, dx ( and dt denote the advection speed, the spatial mesh 

c interval and the time step size, respectively. Let (i) the 

c Courant number cn = a*dt/dx,- and (ii) theta denote the phase 

c angle of a Fourier mode. Then, according to Proposition 1(a), 

c for any pair of cn and theta with |cn| . le. 1/2 and -pi .It. 

c theta .le. pi, the three amplification factors of the a3 scheme 

c are all of unit magnitude, 
c 

c Assuming that |cn| .le. 0.5, this code can be used to evaluate 

c the dispersive errors of the principal amplification factor 

c (per dt) of the a3 scheme and compare them with the 

c corresponding errors of the dual a scheme (see AIAA 2006-4779) . 
c 

c Let lambda = the wavelength of a Fourier mode. Then (i) theta 
c = 2 *pi*dx/ lambda; (ii) the phase angle of the analytical 
c amplification factor is -cn* theta. Let phi denote the phase 

c angle of the principal amplification factor of the dual a scheme 
c or the a3 scheme. Then phi = f(cn, theta) with 

c f(-cn, theta) = f(cn, -theta) = -f(cn, theta) (see Eq.(3.14) in 
C AIAA 2007-4321) . 

c Thus, without any loss of generality, one may assume that 

c 0 .le. cn .le. 0.5 and 0 ,le. theta .le. pi in numerical 

c computations . 

c 

c nz = number of intervals over the domain 0 .le. theta .le. pi. 

c ep = the error bound below which the Newton's iteration is 

c terminated, 

c 

c In the output, (i) cn is the courant number (ii) z is the value 

c of theta, (iii) era is the value of [1-phi/ (-cn*theta) ] for 
c the dual a scheme, (iv) era3 is that for the a3 scheme, and 

c (v) k is number of Newton's interations required for convergence 

c 

c The domain of cn is 0 . le . cn .le. 0.5. 
c 

cn = O.OOldO 

nz = 100 

ep = l.d-14 

pi = 3 . 14 1 5 92 6 53 5 8 97 93 2d0 

dz = pi/df loat (nz) 

z = 0 . dO 
c 

open (unit=8 , f ile= ' fa3 . opt ' ) 
write (8,10) 
write (8,20) cn,nz,ep 
do 300 i = l,nz 
z = z + dz 

a = 2 . d0*cn* *2* ( 1 . dO - dcos(z)) 
b = cn*dsin(z) 
pO = -cn*z 

era = l.dO - datan(-b/dsqrt (1 .dO - b**2})/p0 

p = pO 

k = 0 

100 pn = p - ( (dsin (p) ) **2 - a*(l.d0 + dcos(p)) - b?dsin(p))/ 

* (dsin{2 .d0*p) + a*dsin(p) - b*dcos(p)}« 

k = k + 1 

if (dabs(pn/p - l.dOJ.lt.ep) goto 200 
P = pn 



goto 100 


200 

era3 = l.dO - pn/pO 



write (8,30) z,era,era3,k 

300 

continue 
close (unit=8) 


10 

format (' ************ 

, ' Output for Code fa3.for',' *********** 

20 

format (' cn =',gl4.7, 

' nz =',i6,’ ep =',gl4.7) 

30 

format (' z =',gl4.7,' 

stop 

end 

era =',gl4.7,' era3 =',gl4.7,' k=',i6) 
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Figure 1 . — A surface element on the boundary 
S(V) of an arbitrary space-time volume V. 


j = -5 j = -3 j = -1 j = 1 j = 3 j = 5 



2(a). — The space-time mesh. 



Figure 2. — The SEs and CEs. 


NASA/TM— 2008-215138 


70 



phase error (a-scheme) 



Figure 3. — Phase errors of the a(3) scheme and the dual a scheme. 
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phase error (a3-scheme) 
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